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I.  INTRODUCTION 


A.  DESCRIPTION  OF  THE  PROBLEM 

A  statistical  design  is  a  plan  according  to  which  an  experiment  is  patterned.  It 
provides  the  basis  upon  which  appropriate  statistical  tests  and  inferences  can  be  made 
after  the  experiment  has  been  performed.  The  selection  of  the  experimental  design  to 
be  used  in  a  given  situation  is  extremely  important  because  it  plays  a  predominant  role 
in  the  efficiency  of  the  experiment,  the  precision  with  which  the  objectives  are  met,  and 
the  total  effort  (and  cost)  expended  upon  the  experiment. 

The  concern  of  this  thesis  is  the  development  of  power  algorithms  for  several 
common  tests,  including  the  F-test  and  t-test.  The  experimenter  wants  a  test  to 
achieve  the  correct  decision  with  as  high  probability  as  possible  in  practical  operational 
tests  and  evaluations  given  certain  conditions  hold. 

B.  SCOPE  OF  THE  THESIS 

An  interactive  computer  program  is  developed  which  uses  power  algorithms  for 
several  common  tests.  The  program  can  be  used  by  decision  makers  who  may  not 
have  deep  knowledge  of  statistics  or  may  be  unfamiliar  with  the  details  of  the  design  of 
experiments.  The  output  helps  the  decision  maker  design  experiments  appropriate  for 
testing  hypotheses  with  common  statistical  tests.  The  experiment  designer  can  use  the 
program  to  evaluate  the  effect  upon  power  of  varying  underlying  design  parameters 
such  as  sample  size. 

Because  the  underlying  statistical  noncentral  distributions  used  in  computing 
power  do  not  have  tables,  approximation  methods  are  used  for  computing  the 
noncentral  CDFs.  Approximation  methods  were  developed  which  are  efficient, 
requiring  reduced  computer  CPU  time. 

C.  BACKGROUND 

A  hypothesis  is  a  statement  about  the  values  of  the  parameters  of  a  probability 
distribution.  For  example,  suppose  we  think  that  the  mean  yield  of  a  chemical  process 
is  more  than  94.5  percent.  Hypotheses  to  test  this  statement  might  be  expressed 
formally  as 

Ho  :  p  <  94.5 
Ha  :  p  >  94.5 


9 


The  statement  Ho  :  H  ^  94.5  is  called  the  null  hypothesis,  and  Ha  :  ft  >  94.5  is  called 
the  alternative  hypothesis.  The  value  of  the  mean  specified  in  the  null  hypothesis 
might  typically  be  determined  in  one  of  two  ways.  It  may  be  the  result  of  some  theory 
or  model  regarding  the  process  under  study,  or  it  may  be  the  result  of  contractual 
specifications. 

To  test  a  hypothesis  one  usually  devises  a  procedure  for  taking  a  random  sample, 
computes  an  appropriate  test  statistic,  and  then  rejects  or  declines  to  reject  the  null 
hypothesis  Ho,  depending  on  the  outcome  on  the  test  statistic.  Part  of  this  procedure 
is  the  specification  of  the  set  of  values  for  the  test  statistic  which  lead  to  rejection  of 
Ho.  This  set  of  values  is  called  the  critical  region  or  rejection  region  for  the  test. 

Two  kinds  of  errors  may  be  committed  when  testing  hypotheses.  If  the  null 
hypothesis  is  rejected  when  it  is  true,  then  a  type  I  error  has  occurred.  If  the  null 
hypothesis  is  not  rejected  when  it  is  false,  then  a  type  II  error  has  been  made.  The 
probabilities  of  these  two  errors  are  given  special  symbols: 

«  =  P(type  I  error)  =  P( reject  Ho  |  Ho  is  true) 
p  =  P(type  II  error)=  P(fail  to  reject  Ho  |  Ho  is  false) 

Generally,  the  determination  of  p  requires  the  use  of  a  noncentral  sampling 
distribution  (e.g.,  noncentral  F,  noncentral  /  “  ,  and  noncentral  t)  for  which  tables  are 
not  readily  available.  The  noncentral  sampling  distribution  depends  on  a  parameter 
called  the  noncentrality  parameter.  The  noncentrality  parameter  is  usually  a  measure 
of  the  distance  (in  some  sense)  between  the  values  of  the  parameter  under  the  null  and 
alternate  hypotheses.  Thus,  type  II  error  rates  depend  upon  parameters  of  the 
distribution  of  the  test  statistic  under  the  alternative  hypothesis.  The  power  of  the  test 
for  a  specified  value  of  the  hypothesized  parameter  is  the  probability  the  test  would 
reject  Ho  when  Ho  is  false,  or  1  —  p. 

In  addition,  whatever  the  test  procedure  is,  rejecting  Ho  when  Ho  is  not  true  is 
usually  something  we  want  a  test  to  achieve  with  as  high  a  probability  as  possible. 
Therefore  we  want  the  power  of  the  test,  for  a  given  value  of  non-centrality  parameter, 
to  be  high.  Power  algorithms  are  used  to  assist  in  the  "best"  selection  of  parameters 
for  a  statistical  design.  There  do  not  exist  simple  closed  form  equations  for  computing 
the  required  probabilities  from  the  common  noncentral  sampling  distributions. 
Consequently,  one  must  use  numerical  approximations  to  determine  power.  In  this 
thesis  we  use  various  approximation  methods  to  compute  the  cumulative  probabilities 
for  the  noncentral  F,  noncentral  x  >  and  noncentral  t  distributions. 


10 


II.  NONCENTRAL  CHI-SQUARE, T,F  DISTRIBUTIONS 


A.  NONCENTRAL  CHI-SQUARE  DISTRIBUTION 

If  Wl,  W2,  ....  Wn  arc  independently  distributed  as  N(p-  ,1),  i=  1,2,  ...,n,  then 

/y  '•y 

y  Wiz  has  a  distribution  known  as  a  noncentral  /“  (n,  X  )  where  n  is  the  number  of 
degrees  of  freedom  and  X  =  Y  ^  “  is  the  noncentrality  parameter.  When  fi  |  =  p 

= ...  =  (ii  =  0,  then  X  =  0,  and  the  noncentral  xz  (n,0)  reduces  to  the  usual  central 
X"  (n)  with  n  degrees  of  freedom.  The  cumulative  distribution  function  of  x"  (n,X  )  is, 


-  e 


- 


1 


;Vo  3!  [»(**;) 


*  ; 


(2.1) 

x  x  o 


while  F(x;n,  X)  =  0  for  x  <  0 

It  is  possible  to  express  F(x;n,X)  for  x>0,  in  an  easily  remembered  form  as  a 

q 

weighted  sum  of  central  x^  CDF's  with  weights  equal  to  the  probabilities  of  a  Poisson 
distribution  with  expected  value  X/2.  That  is, 


pU;n,A)  =  £ 

)=o 


(2.2) 


Thus  a  X  (n,  X)  variable  can  be  regarded  as  a  mixture  of  central  x  variables.  This 

2 

interpretation  is  often  useful  in  deriving  the  distribution  of  functions  of  noncentral  x 
random  variables. 

The  probability  density  function  can,  similarly,  be  expressed  as  a  mixture  of 
2 

central  X  probability  density  functions 


where 


(2.3) 
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The  mean  and  variance  of  the  distribution  are 

E(X)=  n  +  2  X  and  Var(X)=  2n+ 8X  [Ref.  l:p.  130], 

B.  NONCENTRAL  F  DISTRIBUTION 

IfYj  and  Y2  are  independent  and  Y]  is  x  j  ^  fa,  ,X  )  and  Y2  is  x  “  (n2  )  then 


V  =  (Y1/nl)/(Y2/n2) 


(2.4) 


is  distributed  as  F'(n1  ,n,  ,X  ),  the  noncentral  F  distribution  with  nt  and  n,  degrees  of 
freedom  and  noncentrality  parameter  X  .  Its  density  function  is  given  by 


00 


.  AVp(i+^tKJ^ 


where 


c  = 


G  nj37 


(2.5) 


when  X.  =0  this  reduces  to  the  density  function  of  the  central  F(nltn2  )  distribution. 
The  mean  and  variance  of  the  distribution  are 


forn>2 


and 


Var(V) 


=  aYi^  r 

7'.*‘(1'ci-)t 


C*r*)(n«r4)  ru-  4- 


(for  n?  >  4) 


When  X  =0,  these  reduce  to  the  mean  and  variance  of  the  central  F(nj,n2) 
distribution.  Derivation  of  (2.5)  is  shown  in  [Ref.  l:p.  189]. 

The  cumulative  distribution  of  V  can  be  expressed  in  terms  of  an  infinite  series  of 
multiples  of  incomplete  beta  functions,  as  follows: 
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is  the  incomplete  beta  function. 


where  I  », 

flit- n,}. 


In  this  thesis,  we  use  Paulson's  Approximation  from  Severo  and  Zelen  [Ref.  2]. 
Their  approximation  is  as  follows: 


Pr(V  <  f  0  )  ~  0)  (x) 


(2.7) 


where 

x  = 


i 


(‘-ft  )(-£■,) 

'-1 

[i- 

An* 

.+X|)  ] 

+f 

■Kl(  n'f°  ) 

>fj 

!4 

and 


<I>  is  the  standard  Normal  CDF. 

The  following  table  shows  some  values  obtained  using  this  approximation, 
together  with  exact  values  of  Pr(V  ^  fo  ).  It  can  be  seen  in  TABLE  1  that  this 
approximation  gives  about  2  decimal  place  accuracy.  Thus,  this  approximation 
provides  adequate  accuracy  for  computing  the  power  of  F-tests  [Ref.  3:p.  202). 

The  program  we  used  in  the  computation  of  noncentral  F  CDF  values  is  shown 
in  Appendix  D. 


C.  NONCENTRAL  T  DISTRIBUTION 

The  ratio 


T  =  (U  +  X)/  V(Y/n)  (2.8) 

where  U  and  Y  are  independent  random  variables  distributed  as  standard  normal 
(N(0,1))  and  x  ^  with  n  degrees  of  freedom  .respectively,  is  said  to  have  a  noncentral  t 
distribution  with  n  degrees  of  freedom  and  noncentrality  parameter  X.  Sometimes  X/ 
(or  even  .5X^),  rather  than  X,  is  termed  the  noncentrality  parameter.  If  X  is  equal  to 
zero,  the  distribution  is  central  t  with  n  degrees  of  freedom. 

The  cumulative  distribution  of  T  is  [Ref.  l:p.  20 lj 
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Pr(T  <  t)  =• 


(2.9) 


TABLE  1 

CDF  OF  A  NONCENTRAL  F  DISTRIBUTION 


nl 

n2 

X 

fo 

approx 

exact 

3 

10 

4 

3.708 

0.7499 

0.745 

4 

6.552 

0.9192 

0.918 

16 

3.708 

0.2023 

0.206 

16 

6.552 

0.5190 

0.517 

3 

20 

4 

3.098 

0.7067 

0.700 

4 

4.938 

0.8894 

0.887 

16 

3.098 

0.1186 

0.126 

16 

4.938 

0.3488 

0.347 

5 

10 

6 

3.326 

0.7337 

0.731 

6 

5.636 

0.9143 

0.914 

24 

3.326 

0.1553 

0.158 

24 

5.636 

0.4629 

0.461 

5 

20 

6 

2.711 

0.6685 

0.664 

6 

4.103 

0.8715 

0.870 

24 

2.711 

0.0643 

0.069 

24 

4.103 

0.2437 

0.245 

8 

10 

9 

3.072 

0.7159 

0.714 

9 

5.057 

0.9087 

0.908 

36 

3.072 

0.1166 

0.119 

36 

5.057 

0.4088 

0.408 

8 

30 

9 

2.266 

0.581 

0.578 

9 

3.173 

0.8157 

0.813 

36 

2.266 

0.0146 

0.017 

36 

3.173 

0.0846 

0.088 

In  this  thesis,  we  use  an  approximation  method  suggested  by  Abamowitz  and 
Stegun  [Ref.  4:p.  949]. 

Pr(T<t)~<D(x)  (2.10) 

where  x  = -  and  <I>  is  the  standard  Normal  CDF. 
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The  following  table  shows  some  values  obtained  using  this  approximation, 
together  with  exact  values  of  Pr(T  ^  to  ).  It  can  be  seen  in  Table  2  that  this 
approximation  gives  about  2  decimal  place  accuracy.  Thus,  this  approximation 
provides  adequate  accuracy  for  computing  the  power  of  t-tests. 

The  program  we  used  in  the  computation  of  noncentral  t  CDF  values  is  shown 
Appendix  E. 


TABLE  2 

CDF  OF  A  NONCENTRAL  T  DISTRIBUTION 


n 

X 

to 

approx 

exact 

4 

0.112 

4 

0.017 

0.01 

7.944 

4 

0.993 

0.99 

2.24 

10 

0.026 

0.01 

18.494 

10 

0.993 

0.99 

24 

5.786 

9.797 

0.012 

0.01 

13.981 

9.797 

0.991 

0.99 

16.118 

24.494 

0.014 

0.01 

33.078 

24.494 

0.992 

0.99 
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III.  POWER  OF  THE  F-TEST 


A.  INTRODUCTION 

Under  appropriate  conditions,  the  best  test  for  testing  equality  of  several  means 
is  the  analysis  of  variance  test.  Analysis  of  variance  has  a  wide  application.  It  is  one 
of  the  most  useful  techniques  in  the  field  of  statistical  inference.  As  in  any 
hypothesis-testing  situation,  the  power  of  the  F  test  is  of  interest  to  the  experimenter. 

In  this  chapter  we  will  discuss  the  power  of  F  tests  and  provide  an  example.  To 
give  an  overall  evaluation  of  the  power  of  F  tests  in  the  analysis  of  variance,  we  may 
use  power  curves.  An  important  use  of  the  power  curve  is  to  guide  the  experimenter  in 
selecting  fhe  sample  size  (number  of  replicates)  so  that  the  design  will  be  sufficiently 
sensitive  to  important  potential  differences  in  the  treatments.  We  will  consider  power 
curves  for  one-way  and  multi-way  analyses  of  variance  (ANOVA's). 

B.  THE  ONE-WAY  CLASSIFICATION  ANALYSIS  OF  VARIANCE 

Suppose  we  have  k  different  levels  of  a  single  factor  that  we  wish  to  compare. 
The  different  levels  of  the  factor  are  often  called  treatments.  The  observed  responses 
from  each  of  k  treatments  is  a  random  sample  on  a  random  variable.  The  data  would 
appear  as  in  Figure  3.1. 


Treatment 


Observation 


1 

Y11 

Y12 

■  Yln 

2 

Y21 

y22 

•  Y2n 

• 

. 

. 

• 

k 

Ykl 

C\J 

>- 

Ykn 

Figure  3.1  Typical  Data  for  One-way  Classification  ANOVA. 
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We  will  find  it  useful  to  describe  the  observations  by  the  linear  statistical  model: 


Y  y  =  H  +  +  £  ij  (3.1) 

where  Y*j  is  the  (ij)^  observation,  p  is  a  parameter  common  to  all  treatments  called 
the  overall  mean,  T-  is  a  parameter  peculiar  to  the  itF  treatment  called  the  itF  treatment 
effect,  and  C-  is  a  random  error  component,  assumed  to  be  I  ID  N  (0,  <sL  ). 

The  indices  used  are: 

•  i  =  the  treatments,  i  =  l,2,...,k 

•  j  =  the  replication  per  treatment,  j  =  l,2,...,n 

The  objective  in  the  ANOVA  is  to  test  appropriate  hypotheses  about  the 
treatment  effects.  The  variance  a~  is  assumed  constant  for  all  levels  of  the  factor. 
This  model  is  called  the  one-way  classification  analysis  of  variance  because  only  one 
factor  is  investigated. 

We  are  interested  in  testing  the  equality  of  the  k  treatment  effects,  so  the 
appropriate  hypotheses  are 

Ho  :  T |  =  =  ...  =  =  T 

Ha  :  Tj  #  t:  for  some  i,j  . 

That  is,  if  the  null  hypothesis  is  true,  then  each  observation  is  made  up  of  the  mean  p 
+  t  plus  a  realization  of  the  random  error  £”. 

The  results  of  the  ANOVA  procedure  is  summarized  in  Table  3. 

We  may  use  the  expected  values  of  the  mean  squares  to  verify  that  Fo  (in  Table 
3)  is  an  appropriate  test  statistic  for  Ho.  From  the  expected  mean  square  we  see  that, 

9 

in  general,  MSE  is  an  unbiased  estimator  of  a  .  Also,  under  the  null  hypothesis,  MSt 

is  an  unbiased  estimator  of  a  .  However,  if  the  null  hypothesis  is  false,  then  the 

expected  value  of  MSt  is  +(n  £  (Tj)")/(k-l).  The  expected  value  of  the  mean 
2 

square  error  is  <r\ 

Therefore,  under  the  alternate  hypothesis  the  expected  value  of  the  numerator  of 
the  test  statistic  (Fo)  is  greater  than  the  expected  value  of  the  denominator  and  we 
would  reject  Ho  with  values  of  the  test  statistic  which  are  too  large.  That  is,  we  would 
reject  Ho  if 

Fo  >  Fa>(k-l)>(N-k) 
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TABLE  3 

THE  ANALYSIS  OF  VARIANCE  FOR  A  ONE-WAY  CLASSIFICATION 


Source  of 
Variation 

Sum  of 
Squares 

Degrees  of 
Freedom 

Mean 

Square 

Fo 

Between 

Treatments 

sst 

k-1 

MSt 

MSt/MSE 

Error(within 

treatments) 

SSE 

k(n-l) 

MSE 

Total 

SST 

kn-1 

where  the  numerator  and  denominator  df  are  k-1  and  k(n-l)  =  N-k,  and  a  is  the  type  I 
error  rate. 

The  power  of  the  test  is: 


1-P  =  P(F0  >  Fa,(k_1)l(N_k)  |  Ho  is  false) 


(3.2) 


To  evaluate  the  P  in  Equation  3.2  we  need  to  know  the  distribution  of  the  test  statistic 
Fq  if  the  null  hypothesis  is  false.  It  can  be  shown  that,  if  Ho  is  false,  the  statistic  FQ 
has  the  noncentral  F  distribution  with  k-1  and  k(n-l)  degrees  of  freedom  and 
noncentrality  parameter  X,  given  by 

«E(Ti)2 

l  =  - 


The  noncentrality  can  be  interpreted  as  the  squared  standardized  distance 

between  the  origin  and  (tj^,  •••  Ak).  The  ratio  Y  (tj  )^  /  <x^  is  called  the  squared 

?  2 

standardized  distance.  If  only  an  estimate  of  G  is  available,  one  may  replace  gl  with 
the  estimate  [Ref.  5:p.  34). 


C.  THE  MULTI-WAY  CLASSIFICATION  ANALYSIS  OF  VARIANCE 

Many  experiments  require  a  study  of  the  effects  of  two  or  more  factors.  It  can 
be  shown  that,  under  certain  conditions,  factorial  arrangements  are  the  most  efficient 
designs  for  this  type  of  analysis. 
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One  of  the  simplest  factorial  experiments  involves  only  two  factors  or  sets  of 
treatments  say  factor  A  and  factor  B.  Suppose  there  are  a  levels  of  factors  A  and  b 
levels  of  factor  B,  and  these  are  arranged  in  a  2-way  factorial  design;  that  is,  each 
replication  of  the  experiment  contains  all  ab  treatment  combinations.  Assume  there 
are  n  replications  of  the  experiment,  and  let  represent  the  observation  taken  under 
the  i^  level  of  factor  A  and  the  j1^1  level  of  factor  B  in  the  k^  replication. 

The  data  can  be  summarized  as  shown  in  Figure  3.2.  The  order  in  which  the  abn 
observations  are  taken  is  selected  at  random. 


Factor  B 


1  1 

1  | 

2  |  . . 

.  |  b 

1  1 

Y111’Y112  1 

Y121 ,Y122  1 

1  Ylbl ,Ylb2 

1  1  1,. 

Y12n  1 

>•••’  Y12n  1 

1 

••  1  «•••>  Ylbn 

1 

1  1 

Y211 ,Y212  1 

Y221 ,Y222  1 

1  Y2bl’Y2b2 

Factor 

1  2  | 

*•••»  Y21n  1 

»• •• ’  Y22n  1 

1 

• •  1  ••• •  >  Y2bn 

1 

A 

...  | 

1 

1 

1 

...1 

1  1 

Yall ,Yal2  1 

Ya21 ’ Ya22  1 

1  Yabl ,Yab2 

1  a  1 

’•  •  •  >  Ya2n  1 

>•••>  Ya2n  1  • 

1 

•  •  1  »•  •  •  ’  Yabn 

1 

Figure  3.2  Typical  Data  notation  for  a  Two-way  Classification.. 

The  observations  may  be  described  by  the  linear  model: 

Yijk  =  B  +  ti  +  pj+  Tp-  +£ijk,  (3.3) 

where 
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•  H  =  the  overall  mean  effect; 

•  t  j  =  the  true  effect  of  i  ^  level  of  factor  A,  i  =  l,2,...,a; 

•  P  j  =  the  true  effect  of  j  ^  level  of  factor  B,  j  =  l,2,...,b; 

•  r  P  ::  =  the  effect  of  the  interaction  between  T  j  and  P  :  ;  and 

•  c  ^  =  a  random  error  component,  assumed  to  be  III)  N  (0,(T  ). 

Both  factors  are  assumed  to  be  fixed.  It  is  usually  assumed  that  the  treatment  effects 
are  defined  as  deviations  from  the  overall  mean,  so  V  T  j  =  0  and  £  p  j  =  0. 
Similarly,  the  interaction  effects  are  fixed  and  usally  defined  so  that  Y,  (  t  p  )  ij  =  0, 
where  the  summation  is  over  either  i  or  j.  Since  there  are  n  replicates  of  the 
experiment,  there  are  a  total  of  abn  observations. 

We  are  interested  in  testing  various  hypothesis  about  the  parameters  in  equation 
3.3  .  An  appropriate  hypothesis  testing  procedure  would  again  be  analysis  of  variance. 
More  specifically,  as  we  are  considering  two  controllable  sources  of  variation  (A  and 
B),  the  procedure  is  called  the  two-way  classification  analysis  variance. 

In  order  to  test  the  hypothesis  Ho  ;  t  j  =  0  for  i=l,2,...,a  (no  row  factor 
effects),  Ho  ;  p  j  =  0  for  j=  l,2,...,b  (no  column  factor  effects),  and  Ho  ;  (  X  p  )  -  =  0 
(no  interaction  effects),  we  can  express  the  total  sum  of  square  as: 

SST  =  SSA  +  SSB  +  SSAB  +  SSE.  (3.4) 

Here,  SSA  is  a  sum  of  squares  due  to  "rows"  or  factor  A,  SSB  is  a  sum  of  squares  due 
to  "columns"  or  factor  B,  SSAB  is  a  sum  of  squares  due  to  the  interaction  between  A 
and  B,  and  SSE  is  a  sum  of  squares  due  to  error.  The  degrees  of  freedom  associated 
with  each  sum  of  squares  are  shown  in  Figure  3.3. 

If  we  assume  £  ^  are  IID  X(0,<7"  )  and  apply  Cochran's  theorem  (Theorem  3-1) 
under  the  null  hypothesis  of  no  effects,  each  sum  of  squares  on  the  right-hand  side  of 
Table  4  when  divided  by  <r  is  distributed  as  x  with  the  indicated  number  of  degrees 
of  freedom,  and  these  statistics  are  independent. 

Theorem  3-1  (COCHRAN).  Let  Z  j  be  IID  N(0,1)  for  i=  l,2,...,v  and  suppose 
Yz  j  ^  =  Qj  +  q9  4-  q3  -k..  +  Qs  where  s  ^  v,  and  the  quadratic  form  Qj  has 
v-  degrees  of  freedom  (i=  l,2,...,s).  Then  Qj  Q2  Qs  are  independent 
chi-square  random  variable  with  Vj  V2  ..,  vs  degrees  of  freedom,  respectively,  if 
and  only  if 

v  =  yl  +  v2  +  ...+  vs. 
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Effect 

Degrees  of  freedom 

A 

a3! 

B 

b-1 

AB  interaction 

(a-l)(b-l) 

Error 

ab(n-l) 

Total 

abn-1 

Figure  3.3  Table  of  degrees  of  freedom  with  sums  of  squares.. 

Assuming  that  factors  A  and  B  are  fixed,  the  expected  values  of  the 
mean  squares  are: 

E(MSA)  =  <r2  +(  bn  51  t  j  2  )/(a-l) ;  (3.5) 

E(MSB)  =  <T2  +(bn£pi2)/(b-l);  (3.6) 


E(MSAB)  =  <r2  +  (n  £  £  (t  P  )  -j  2  )/  (a-l)(b-l) ;  and  (3.7) 

E(MSE)  =  <T2  .  (3.8) 

Therefore,  to  test  the  hypothesis  Ho  :  t  •  =  0  ;  i  =  l,2,...,a  (no  row  factor  effects),  and 
Ho  ;  P  j  =  0  ;  j  =  1,2,. ..,b  (  no  column  factor  effects),  and  Ho  :  (  r  p  )  jj  =  0  (no 
interaction  effects),  we  would  divide  the  corresponding  mean  square  by  the  mean 
square  error.  Under  the  null  hypothesis  of  no  effect,  this  ratio  will  follow  an  F 
distribution  with  appropriate  numerator  degrees  of  freedom  and  ab(n-l)  denominator 
degrees  of  freedom,  and  the  critical  region  will  be  located  in  the  upper  tail.  The  test 
procedure  is  usually  summarized  in  an  analysis  of  variance  table,  such  as  shown  in 
Table  4. 
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TABLE  4 

THE  ANALYSIS  OF  VARIANCE  FOR  A  TWO-WAY  CLASSIFICATION 


Source  of 
Variation 

Sum  of 
Squares 

Degrees  of 
Freedom 

Mean 

Square 

Fo 

A  treatments 

SSA 

a-1 

MSA 

MSA/MSE 

B  treatments 

SSB 

b-l 

MSB 

MSB/MSE 

Interaction 

SSAB 

(a-l)(b-l) 

MSAB 

MSAB/MSE 

Error 

SSE 

ab(n-l) 

MSE 

Total 

SST 

abn-1 

To  the  compute  the  power  of  tests  in  two-way  ANOVAs.  the  procedure  is  the 
same  as  in  the  one-way  case.  Relationships  among  X,  the  numerator  degrees  of 
freedom  and  the  denominator  degrees  of  freedom  are  shown  in  Table  5  . 

In  a  similar  way  one  can  expand  to  the  general  multi-way  AN'OVA  procedure 
[Ref.  5:p.  124], 

D.  THE  ALGORITHMS  AND  FLOWCHART 

A  program  for  computing  a  power  table  and  power  curve  in  general  ANOVA's 
(fixed  model)  is  shown  in  Appendix  F.  The  power  of  one-way  ANOVA's  are  solved 
with  these  programs,  using  the  following  the  sequence: 

1)  Determine  what  variable  to  include,  such  as  sample  size  vs  power,  number  of 
treatments  vs  power,  a-level  vs  power,  or  noncentrality  vs  power. 

2)  Determine  the  values  of  the  maximum,  minimum,  and  increment  of  the  variable 
chosen. 

3)  Given  input  data,  compute  the  critical  value,  using  the  inverse  of  the  central 
F-distribution  (Appendix  A). 

4)  Compute  the  power  value  in  each  case;  that  is,  the  CDF  of  the  noncentral 
F-distribution  (Chapter  II). 

5)  Print  the  power  table  or  the  power  curve. 
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TABLE  5 

NONCENTRALITY  PARAMETERS  FOR  POWER  IN  A  TWO-WAY  ANOVA 


Factor  X 

Numerator 

DOF 

Denomi nator 

DOF 

A 

a-1 

ab(n-l) 

B 

b-1 

ab(n-l) 

AB 

Pjj2 

"A 

(a-l)(b-l) 

ab(n-l) 

Multi-way  ANOVAs  have  the  same  algorithms  but  multi-way  ANOVA  may 
include  tests  of  interaction  efFects.  The  program  considers  only  up  to  three-way 
interaction  effects.  We  can  explain  how  to  compute  degrees  of  freedom  of  error  term 
in  m-way  ANOVA's,  assuming  the  balanced  case,  as  follows:  Total  degrees  of  freedom 
is  DOF(Total)  =  n  £  k(i)  -1,  where  k(i)  is  the  number  of  levels  of  the  i*^  factor  [Ref.  9]. 

(1)  If  the  model  has  only  main  effects  without  any  interaction  effects, 

DOFl(Error)=  DOF(Total)  -  £  (k(i)-l). 

(2)  If  the  model  has  several  factors  and  only  2-way  interaction  effects, 

DOF2(Error)=  DOFl(Error)  -  £  £  (k(i)-l)(k(j)-l). 

(3)  IF  the  model  has  several  factors  and  only  3- way  interaction  effects, 

DOF3(Error)=  DOFl(Error)  -  V  £  £  (k(i)-l)(k(j)-l)(k(k)-l). 

(4)  If  the  model  has  several  factors  and  2-way  and  3-way  interaction 

effects, 

DOF4(Error)  =  DOF3(Error)  -  £  £  (k(i)-l)(k(j)-l). 

If  the  model  has  more  than  3-way  interaction  effects,  then  the  user  must  modify 
the  degrees  of  freedom  of  the  error  terms  accordingly.  A  flowchart  is  shown  in  Figure 
3.4. 
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Figure  3.4  System  flowchart  for  F-test  power. 
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E.  EXAMPLES  OF  F-TESTS 


1.  Example  of  One-way  ANOVA  test  (number  of  replicates  vs  power) 

a.  Scenario 

A  manufacturer  suspects  that  the  batches  of  raw  material  furnished  by  his 
supplier  differ  significantly  in  calcium  content.  There  are  a  large  number  of  batches 
currently  in  the  warehouse.  Five  of  those  are  randomly  selected  for  study.  A  chemist 
wants  to  know  the  appropriate  sample  size  per  each  batch,  in  order  to  test 
Ho  1  Tj  =  =  •••  = 

against 


Ha  :  Tj  #  Tj  for  some  i,j 


The  chemist  would  like  to  know  how  many  replicates  to  run  if  it  is  important  to  reject 
Ho  with  probability  at  least  0.9  when  the  standardized  distance  square 
is  2  and  (l  =  .05.  Thus  he  would  like  to  know  what  the  power  of  the  F-test  is  for  a 
range  of  possible  replicates.  He  decides  to  check  replicates  from  2  to  12  in  increments 
of  1. 

b.  Inputs 

1)  Select  one-way  ANOVA  test  (the  number  of  replicates  vs  power). 

2)  The  number  of  treatment  =  5. 

3)  a  -level  =  .05 

4)  Standardized  distance  square  =  2. 

5)  Mrximum  replicates  is  12,  minimum  replicates  is  2,  increment  is  1. 

c.  Start  program 

Do  you  want  to  analyze  one  way  ANOVA(y/n)? 

y 

Do  you  want  to  plot  n  (#  of  observation  per  treatment  vs  power)(y/n)? 

y 

Number  of  k(  #  of  treatment)? 

? 

5 

Alpha-  level? 

? 

.05 


Standardized  squared  distance  value  ? 

? 
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2. 

Maximum  n  value  (the  maximum  value  on  X-axis)? 

? 

12 

Minimum  n  value  (the  minimum  value  on  X-axis)? 
(condition  :N  must  be  more  than  2) 

? 

2 

Increment  n  value  ? 

? 


d.  Output 

The  output  is  presented  in  Table  6 

TABLE  6 

OUTPUT  OF  THE  ONE-WAY  ANOVA  EXAMPLE  (REPLICATES  VS  POWER) 


DOF1 

DOF2 

a 

F-INVERSE 

X 

POWER 

4 

5 

0.05 

5.29087 

4.00 

0.15373 

4 

10 

0.05 

3.52496 

6.00 

0.30024 

4 

15 

0.05 

3.10397 

8.00 

0.44935 

4 

20 

0.05 

2.91676 

10.00 

0.58677 

4 

25 

0.05 

2.81108 

12.00 

0.70332 

4 

30 

0.05 

2.74323 

14.00 

0.79562 

4 

35 

0.05 

2.69600 

16.00 

0.86449 

4 

40 

0.05 

2.66122 

18.00 

0.91328 

4 

45 

0.05 

2.63455 

20.00 

0.94631 

4 

50 

0.05 

2.61345 

22.00 

0.96776 

4 

55 

0.05 

2.59634 

24.00 

0.98118 

#  of  replicate  POWER 


NN  = 

2. 

POWER  = 

0.15373 

NN  = 

3. 

POWER  = 

0.30024 

NN  = 

4. 

POWER  = 

0.44935 

NN  = 

5. 

POWER  = 

0.58677 

NN  = 

6. 

POWER  = 

0.70332 
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NN  = 

7. 

POWER  = 

0.79562 

NN  = 

8. 

POWER  = 

0.86449 

NN  = 

9. 

POWER  = 

0.91328 

NN  = 

10. 

POWER  = 

0.94631 

NN  = 

11. 

POWER  = 

0.96776 

NN  = 

12. 

POWER  = 

0.98118 

0.98118  +-+ - + - + - + - + - + — * - *-+ 

POWER  I  *  I 


0. 56746 


+ 


* 


0. 15373 

+  -* - 

1 

1 

1 

+ 

1 

1 

1 

1 

1 

1 

1 

1 

1 

+ 

1 

1 

1 

1 

1 

1 

1 

1 

1 

+ 

1 

+ 

2.000 

7.000 

#  OF  REPLICATES  12.00 

X-SCALE:  0.125E  +  00  UNITS 

Y-SCALE:  "\"=  0.138E-O1  UNITS 
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2.  Example  of  one-way  ANOVA  test  (noncentrality  vs  power) 

a.  Scenario 

Five  brands  of  batteries  arc  under  study.  It  is  suspected  that  the  life  (in 
weeks)  of  the  five  brands  is  different.  Five  batteries  of  each  brand  are  tested.  A 
manufacturer  wants  to  know  the  power  as  a  function  of  the  standardized  squared 
distances.  That  is,  would  like  to  know  what  the  power  of  the  F-test  is  for  a  range  of 
possible  standardized  distances  square.  He  decides  to  check  standardized  squared 
distances  from  1  to  5  in  .2  increment,  where  n  =  5  and  a  =  .05. 

b.  Inputs 

1)  Select  one-way  ANOVA  test  (the  standardized  distance  square  vs  power). 

2)  Number  of  treatment  =  5. 

3)  Number  of  replicates  =  5. 

4)  a  -level  =  .05. 

5)  Maximum  standardized  distance  square  is  5  and  minimum  standardized  distance 
square  is  1  and  the  increment  is  .2. 

c.  Start  program 
Start  program: 

Do  you  want  to  analyze  one-way  ANOVA(y/n)? 

y 

Do  you  want  to  plot  n(#  of  observation  per  treatment)  vs  power(y/n)? 

n 

Do  you  want  to  plot  alpha-level  vs  power(y/n)? 
n 

Do  you  want  to  plot  noncentrality  vs  power(y/n)? 

y 

#  of  observations  per  treatment  (n=  )  ? 

? 

5 

#  of  treatment  (k=  )? 

7 

5 

Alpha-level  ? 

? 

.05 

Maximum  standardized  squared  distance  value  range  ? 
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? 

5 

Minimum  standardized  squared  distance  value  range  ? 

? 

1. 

Increment  standardized  squared  distance  value  range  ? 


.2 

d.  Output 

The  screen  output  (  Table  7  )  is  as  follows: 


TABLE  7 

OUTPUT  OF  THE  ONE-WAY  ANOVA  EXAMPLE  (NONCENTRALITY  VS 

POWER) 


DOF1 

DOF2 

a 

F-INVERSE 

POWER 

4 

20 

0.05 

2.91676 

5.00 

0.30382 

4 

20 

0.05 

2.91676 

6.00 

0.36314 

4 

20 

0.05 

2.91676 

7.00 

0.42203 

4 

20 

0.05 

2.91676 

8.00 

0.47946 

4 

20 

0.05 

2.91676 

9.00 

0.53461 

4 

20 

0.05 

2.91676 

10.00 

0.58677 

4 

20 

0.05 

2.91676 

11.00 

0.63550 

4 

20 

0.05 

2.91676 

12.00 

0.68050 

4 

20 

0.05 

2.91676 

13.00 

0.72163 

4 

20 

0.05 

2.91676 

14.00 

0.75885 

4 

20 

0.05 

2.91676 

15.00 

0.79223 

4 

20 

0.05 

2.91676 

16.00 

0.82192 

4 

20 

0.05 

2.91676 

17.00 

0.84812 

4 

20 

0.05 

2.91676 

18.00 

0.87108 

4 

20 

0.05 

2.91676 

19.00 

0.89107 

4 

20 

0.05 

2.91676 

20.00 

0.90835 

4 

20 

0.05 

2.91676 

21.00 

0.92321 

4 

20 

0.05 

2.91676 

22.00 

0.93591 

4 

20 

0.05 

2.91676 

23.00 

0.94672 

4 

20 

0.05 

2.91676 

24.00 

0.95586 

4 

20 

0.05 

2.91676 

25.00 

0.96356 
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POWER 


NONCENTRALITY 


NN=  5.0 

POWER- 

0.30382 

NN  =  6.0 

POWER- 

0.36314 

NN=  7.0 

POWER- 

0.42203 

NN  =  8.0 

POWER- 

0.47946 

NN=  9.0 

POWER- 

0.53461 

NN  =  10.0 

POWER- 

0.58677 

NN=  11.0 

POWER- 

0.63550 

o 

II 

£ 

POWER- 

0.68050 

NN=  13.0 

POWER- 

0.72163 

NN-  14.0 

POWER- 

0.75885 

NN-  15.0 

POWER- 

0.79223 

NN-  16.0 

POWER- 

0.82192 

NN=  17.0 

POWER- 

0.84812 

NN-  18.0 

POWER- 

0.87108 

NN-  19.0 

POWER- 

0.89107 

NN—  20.0 

POWER- 

0.90835 

NN-  21.0 

POWER- 

0.92321 

NN=  22.0 

POWER- 

0.93591 

NN=  23.0 

POWER- 

0.94672 

NN=  24.0 

POWER- 

0.95586 

NN=  25.0 

POWER- 

0.96356 

K-~ 

if 
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0.96356 

POWER 


+-+ 


+ - + - + - + - + - * — *-+ 


* 


* 


0.63369 


+ 


* 


0.30382 

+  -* - +  - 

5.000 

15.00 

NONCENTRALITY 

X-SCALE:  ' 

0.250E  +  00  UNITS 

Y-SCALE:  " 

|"=  0.1 10E-01  UNITS 
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IV.  POWER  OF  THE  T-TEST 


A.  INTRODUCTION 

In  testing  the  hypothesis  that  the  mean  of  a  normal  distribution  is  equal  to  ftQ 
when  the  standard  deviation  is  unknown,  the  t  statistic  may  be  used.  This  test-statistic 
is  a  function  of  the  sample  mean,  the  sample  standard  deviation,  and  the  sample  size. 

The  t-test  has  a  very  wide  application  in  research  areas.  For  example,  we  may 
wish  to  know  whether  product  A  is  better  than  product  B;  or  whether  the  outputs  of 
machines  C  and  D  form  a  homogeneous  mass  of  product;  and  so  forth.  To  answer 
such  questions  we  can  employ  the  t-test  which  is  one  of  the  most  useful  techniques  in 
the  field  of  statistical  inference. 

In  such  hypothesis-testing  situations  the  power  of  t-tests  is  usually  of  interest  to 
the  experimenter.  We  may  use  the  power  curve  to  evaluate  the  power  of  the  t-test. 
The  required  sample  size  might  be  determined  by  referring  to  a  set  of  power  curves. 
We  will  consider  the  power  curves  for  the  one-sample  and  two-sample  cases.  [Ref.  10] 

B.  ONE-SAMPLE  T-TEST 

We  assume  that  a  random  sample,  XI,  X2,  X3,  ...,  Xn  of  size  n  is  taken  from  a 
normal  population  with  mean  ft  and  standard  deviation  <J.  Suppose  we  wish  to  test  the 
hypothesis  that  ft  ^  ftQ  against  the  alternative  hypothesis  that  ft  >  ft0  .  The 
procedure  is  to  reject  the  hypothesis  Ho  :  ft  <  ftQ  if  ((X-  ft0)Vn)/S>  ta  ,  and  to 
accept  Ho  if  otherwise  where  X  and  S  are  the  sample  mean  and  sample  standard 
deviation  and  where  t  a  is  the  1  -  a1*1  percentage  point  of  the  Student  t-distribution 
with  n-1  degrees  of  freedom,  i.e., 

P(  Student  t  >  ta  )  =  a  (4.1) 

With  critical  region  (  ta  ,  oo  ),  the  hypothesis  ft  =  ft0  is  rejected  with  probability 
a  when  Ho  is  true.  If  the  mean  of  the  normal  distribution  is  ft,  >  ftQ  ,  then  the 
power  of  the  above  test  is 

P(  noncentral  t  >  t  a  |  6  =  ((  ft  ,  -  ftQ)  Jn  )/  <r)  =  1-  P  (4.2) 
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where  8  is  the  noncentrality  parameter  for  the  noncentral  t-distribution  with  n-1 
degrees  of  freedom  [Ref.  6:p  .320]. 

C.  TWO-SAMPLE  T-TEST 

Suppose  we  are  given  two  random  samples  from  two  normal  populations  XI,  X2, 
...,  Xn  and  Yl,  Y2,  ....  Ym,  where  the  X's  and  Y's  are  independent  with  mean  values 
Mx  =  Hy  and  common  unknown  variances  <rx  =  <Ty  =  <r“.  Suppose  we  wish  to 
test  the  null  hypothesis  Ho  :  px  ^  My  against  the  alternative  hypothesis  Ha  that  mx 

>  pv.  The  procedure  is  to  reject  the  hypothesis  Ho  if  (X-  Y)  /  J  (  Sp  ^  (1/n  +  I/m)) 

>  t  a  where  t  a  is  a  critical  value  of  the  Student  t-distribution  based  on  m+n—  2 
degrees  of  freedom.  In  this  case  take  8  =  (  Mx  ~  My  )  /  V  ( I/n  +  1/m)  so  the 
probability  of  rejecting  the  hypothesis  when  Mx  >  My  is  equal  to 

P(  noncentral  t>  ta|8)  =  I-P  (4.3) 

where  8  is  the  noncentrality  parameter  for  the  noncentral  t-distribution  with  m+n-2 
degrees  of  freedom  [Ref.  6:p  .321]. 

D.  TWO-SIDED  T-TEST 

For  tests  where  the  alternative  hypothesis  specifies  that  M  ^  Mo  *n  the  test  of 
section  B  or  mx  ^  M  y  in  the  test  of  section  C,  where  the  direction  can  be  either  up  or 
down,  i.e,  m  >  Mo  or  M  <  Mo  and  Mx  >  My  or  Mx  <  My>  we  htive  wrhat  is  known  as  a 
tw'o-sided  test.  In  this  case  the  rejection  procedure  for  the  test  of  section  B  specifies 
that  Ho  is  rejected  if  either  ((X-  M0)Vn)/S  >  t  ^2  or  ((X-  M0  )  V  n  )  /  S  <  -  t  ^2 
and  the  power  of  the  test  is 

P(  noncentral  t  >  t  ^2  \  5  )  +  P(  noncentral  t  <  -t  ^2  |  8  )  =  1  -  P  (4.4) 

where  8  =  (  Mj  -Mq  )  V  n  /  cr  [Ref.  6:p  .323]. 


E.  THE  ALGORITHMS  AND  FLOWCHART 

Listings  of  programs  for  generating  the  power  table  and  curve  for  t-tests  are 
given  in  the  Appendix  E. 

The  following  steps  are  used  to  compute  the  power: 

1)  Determining  whether  one-sample  or  two-sample  t-tests  should  be  used. 
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2)  After  determining  the  test,  give  the  input  data,  (n,  a  ,  fiQ  ,  <r  ) 

3)  Determine  whether  the  one-sided  or  two-sided  test  should  be  used. 

4)  Compute  the  critical  region  (inverse  t-values).  (The  computation  method  is 
provided  in  Appendix  B) 

5)  Compute  the  power  corresponding  to  the  chosen  user  option.  (The 
computation  method  is  discussed  in  Chapter  II) 

6)  Plot  the  power  corresponding  to  each  above  option. 

A  flowchart  is  shown  in  Figure  4.1a 


F.  EXAMPLES  OF  T-TESTS 

1.  Example  of  one  sample  one-sided  t-test 

a.  Scenario 

The  supervisor  in  an  electronic  instrument  plant  is  concerned  with  the 
repair  time  of  an  electronic  equipment.  Standards  specify  that  the  mean  repair  time 
must  be  at  most  40  hours.  For  power  computation,  it  is  assumed  that  the  variance  of 
repairing  time  is  4.0.  The  appropriate  hypothesis  are 
Ho  :  fi  >  40, 

Fla  :  fi  <  40. 

The  experimenter  decides  to  use  a  sample  of  n  =  25  observations  and  a  =  .05.  He  would 
like  to  know  what  the  power  of  the  t-test  is  for  a  range  of  possible  population  means. 
He  decides  to  check  means  from  38  to  40  in  .1  increments. 

b.  Inputs 

INPUTS:  step  1:  one  sample  one-sided  t-test 

step  2:  fi0  =  40 
step  3:  c  =  2.0 
step  4:  a  =  .05 
step  5:  sample  size  =  25 

step  6:  maximum  fij  =  40.,  minimum  fij  =  38,  increment=  .1 

c.  Start  program 
START  PROGRAM: 


Do  you  want  to  test  one  sample  t-test(y/n)? 

y 

Sample  size  ? 

? 


25 

Alpha-level  ? 
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Figure  4.1a  System  flowchart  for  t-test  power. 


35 


Two  saraole  r-test 


Input  :  saziple  size  j 
alpha- "level  j 


Stop 


Figure  4.1b  System  flowchart  for  t-test  power(cont'd 
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JL 


:  newer  table 


A L 


Output  :  plot  the  cur^’e 


Analyze  again  ? 


N 


c 


Stop 


Figure  4.1c  System  flowchart  for  t-test  power(cont'd.) 
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38 


7 


.05 

Mu-zero  value  ? 

? 

40 

Standard  deviation  value  ? 

? 

2 

Do  you  want  to  use  one-sided  test(y/n)? 

y 

Do  you  want  to  test  Ho  :  mu  <  muO  vs  Ha  :  mu  >  muO  (y/n)? 
n 

Do  you  want  to  test  Ho  :  mu  ^  muO  vs  Ha  :  mu  <  muO  (y/n)? 

y 

Maximum  mul  value  (the  maximum  value  on  X-axis)? 

? 

40 

Minimum  mul  value  (the  minimum  value  on  X-axis)? 

(condition  :  mul  must  be  less  than  muO  ) 

7 

38 

Increment  mul  value  ? 

7 

.1 

d.  Output 

The  screen  output  is  as  follows: 
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»*1 

power 

NN=  38.000 

POWER  = 

0.99933 

NN=  38.100 

POWER  = 

0.99849 

NN=  38.200 

POWER  = 

0.99677 

NN=  38.300 

POWER  = 

0.99346 

NN=  38.400 

POWER  = 

0.98742 

N\  =  38.500 

POWER  = 

0.97706 

NN=  38.600 

POWER  = 

0.96028 

NN=  38.700 

POWER  = 

0.93465 

NN=  38.800 

POWER  = 

0.89772 

NN=  38.900 

POWER  = 

0.84755 

NN=  39.000 

POWER  = 

0.78326 

NN=  39.100 

POWER  = 

0.70559 

NN=  39.200 

POWER  = 

0.61711 

NN=  39.300 

POWER  = 

0.52203 

NN=  39.400 

POWER  = 

0.42564 

NN=  39.500 

POWER  = 

0.33355 

NN=  39.600 

POWER  = 

0.25053 

NN=  39.700 

POWER  = 

0.17996 

NN=  39.800 

POWER  = 

0.12338 

NN=  39.900 

POWER  = 

0.08062 

40 


0. 99933 
power 


+-* — * — * — *+ — * — * 


+ 


+-+ 


-+ - + 

*  * 

* 

* 


>+ 


* 


* 


* 


0. 53997  + 


0.  80616E-01  +-+ - + - + - +— 

38.00  38.95 

X-SCALE:  0.237E-01  UNITS 

Y-SCALE:  T=  0.153E-01  UNITS 


* 


+ 


* 


* 


* 


*  I 

*  | 

,+ - + - *-+ 

39. 90 


2.  Example  of  One  sample  two-sided  t-test 
a .  Scenario 

The  supervisor  in  a  rocket  propellant  plant  is  concerned  with  the  burning 
rate  of  a  rocket  propellent.  Standards  specify  that  the  mean  burning  rate  must  be  40 
inches  per  second.  For  power  computation,  it  is  assumed  that  <s  =  4.0.  The 
appropriate  hypotheses  are 

Ho  :  H  =  40 
Ha  :  n  ¥*  40 
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The  engineer  decides  to  use  a  sample  of  n=  25  observations  and  a  =  .05.  He  would  like 
to  know  the  power  of  the  t-test  is  for  a  range  of  possible  population  means.  He 
decides  to  check  means  from  38  to  42  in  .1  increments. 

b.  Inputs 


INPUTS: 


step  1:  one  sample  two-sided  t-test 

step  2:  =  40 

step  3:  <T  =  2.0 

step  4:  a  =  .05 

step  5:  sample  size  =  25 

step  6:  maximum  =42,  minimum  jij  =38,  increment  =  .1 


c.  Start  program 
START  PROGRAM: 

Do  you  want  to  test  one  sample  t-test(y/n)? 


y 

Sample  size  ? 

? 


25 

Alpha-level  ? 

7 

.05 

Mu-zero  value  ? 

7 


40 

Standard  deviation  value  ? 

? 


2 

Do  you  want  to  use  one-sided  test(y/n)? 
n 

Do  you  want  to  use  two-sided  test(y/n)? 

y 

Maximum  mul  value  (the  maximum  value  on  X-axis)? 

7 


42 

Minimum  mul  value  (the  minimum  value  on  X-axis)? 
(condition  :  mul  must  be  less  than  mu- zero) 


42 


? 


38 

Increment  mul  value  ? 

? 

.1 

d.  Output 

The  screen  output  is  as  follows: 


Mj  POWER 


NN=  38.000 

POWER  = 

0.99770 

NN=  38.100 

POWER  = 

0.99525 

NN=  38.200 

POWER  = 

0.99073 

NN=  38.300 

POWER  = 

0.98279 

NN=  38.400 

POWER  = 

0.96965 

NN=  38.500 

POWER  = 

0.94910 

NN=  38.600 

POWER  = 

0.91875 

NN=  38.700 

POWER  = 

0.87639 

NN=  38.800 

POWER  = 

0.82057 

NN=  38.900 

POWER  = 

0.75109 

NN  =  39.000 

POWER  = 

0.66944 

NN=  39.100 

POWER  = 

0.57882 

NN=  39.200 

POWER  = 

0.48379 

NN=  39.300 

POWER  = 

0.38975 

NN=  39.400 

POWER  = 

0.30192 

NN=  39.500 

POWER  = 

0.22459 

NN=  39.600 

POWER  = 

0.16066 

NN=  39.700 

POWER  = 

0.11146 

NN-  39.800 

POWER  = 

0.07708 

NN=  39.900 

POWER  = 

0.05692 

NN=  40.000 

POWER  = 

0.05028 

NN=  40.100 

POWER  = 

0.05687 

NN=  40.200 

POWER  = 

0.07698 

NN=  40.300 

POWER  = 

0.11131 

NN=  40.400 

POWER  = 

0.16045 
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NN=  40.500 

POWER  = 

0.22433 

NN=  40.600 

POWER  = 

0.30162 

NN=  40.700 

POWER  = 

0.38942 

NN=  40.800 

POWER  = 

0.48344 

NN=  40.900 

POWER  = 

0.57848 

NN=  41.000 

POWER  = 

0.66912 

NN=  41.100 

POWER  = 

0.75081 

NN=  41.200 

POWER  = 

0.82034 

NN=  41.300 

POWER  = 

0.87621 

NN=  41.400 

POWER  = 

0.91862 

NN=  41.500 

POWER  = 

0.94901 

NN=  41.600 

POWER  = 

0.96959 

NN=  41.700 

POWER  = 

0.98275 

NN=  41.800 

POWER  = 

0.99070 

NN  =  41.900 

POWER  = 

0.99524 
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0.99770 

power 


0. 52398 


0.  50275E 


-  + _ **_* _ |_ 

** 

* 

* 


+ 


+ 


*  9c  9c  9c 

01  +-+ - + - + - *+** - + - + - +-+ 

38.00  39.95  \x{  41.90 

X-SCALE:  =  0.487E-01  UNITS 

Y-SCALE:  T=  0.158E-01  UNITS 
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V.  DESCRIPTION  OF  THE  POWER  PROGRAM 


A.  PROCEDURE  OVERVIEW 

The  interactive  program  included  in  this  thesis  is  written  in  FORTRAN  77.  It 
was  written  for  uce  on  an  IBM370  from  IBM3278  terminal.  It  is  an  interactive 
program.  The  command  'NONCENT'  is  all  that  is  required  to  start  this  program.  The 
executive  file  name  initializes  the  virtual  machine  environment  and  asks  the  user 
questions  about  the  choice  of  program  and  compilation  requirements.  It  then  activates 
the  selected  program.  Both  ANOVA  and  TTEST  are  interactive  programs  which  will 
be  discussed  in  detail  in  later  chapters.  The  output  from  the  programs  is  presented  on 
the  terminal  screen. 

B.  INSTRUCTIONS  FOR  PROGRAM  ACCESS 

To  start  the  program,  make  sure  you  have  loaded  (NONCEN,  ANOVA,  and 
TTEST)  on  your  disk.  In  CMS  (operating  system  mode),  type  'NONCEN'  and  you  will 
see: 

Please  provide  the  FILENAME  for  your  VS  FORTRAN  program. 

Now  type  the  program  name  you  want,  for  example  'TTEST',  the  response  is: 

Do  you  need  to  compile  your  program  ?  (y/n) 

If  you  want  to  run,  type  'Y'  and  your  program  will  be  loaded. 

A  detailed  description  of  the  ANOVA  (Analysis  of  Variance  Test  )  is  given  in 
Chapter  III  and  of  the  t-test  in  Chapter  IV. 

Following  the  screen  output  from  the  selected  test,  the  user  will  be  asked  some 
questions  about  the  output: 

Do  you  wish  to  BROWSE  your  output?  (Y) 
n 

Print  your  output  file?  (Y) 

n 

Do  you  wish  to  XEDIT  the  program  file?  (Y/N) 

n 

Do  you  wish  to  run  the  program  again?  (Y) 

n 

Then  return  to  CMS  mode. 
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VI.  SUMMARY  AND  CONCLUSIONS 


Power  considerations  are  useful  in  the  design  and  assessment  of  statistical  tests. 
The  calculation  of  power  usually  involves  noncentral  distributions  for  which  tables  of 
probability  are  not  available.  A  search  was  made  for  algorithms  approximating  the 
noncentral  t,  F  and  x  distributions.  Algorithms  giving  sufficient  accuracy  and 
making  efficient  use  of  computer  resources  have  been  implemented  in  this  thesis. 
Listings  of  the  FORTRAN  code  for  these  implementations  are  included. 

An  interactive  program  to  compute  and  display  power  curves  for  several  t-test 
and  F-test  situations  has  been  developed.  This  program  is  user  friendly  and  is 
described  in  this  thesis.  A  listing  of  the  program  is  provided.  This  program  should  be 
useful  to  researchers,  experiment  designers  and  statisticians. 
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APPENDIX  A 

AN  INVERSE  CENTRAL  F  APPROXIMATION 


Let  P  be  the  central  F  complementary  CDF,  and  let  Q  be  the 
complementary  standard  normal  CDF  [Ref.  4: p.  947]. 

Define  yp  by  Q(yp)=p,  and  Fp  by  P(Fp  |  v2  ,v2  )=p. 

Then  Fp  ~  exp(2w) 
where  w  =  wl-(w2)(w3) 

wl  =  (yp  V  (h+X))/h 

w2  =  (l/(v2  -D-l/Cvj  -1)) 

w3  =  (X+(5/6)+2/(3h)) 

h  =  2/((l/(v1  -l)+l/(v2  -1)) 

X=(yp  2  —3 )/ 6 

The  program  we  used  in  the  computation  of  critical  values  of  F-tests  is  shown 
below. 


*  APPROXIMATION  TO  THE  F-  INVERSE  DISTRIBUTION.  * 

*  CALCULATE  F-INVERSE  GIVEN  DOF1 >DOF2 , ALPHA .  * 

*  Q(  F-ALPHA  | DOF1 >DOF 2  )=  ALPHA.  * 

*  * 


REAL  ALPHA >NUMER,DENUM, LAM >XP,YP,FINVER>A,B 
DATA  NF1,NF2, ALPHA/2, 4,. 01/ 

DATA  CO, Cl, C2/2. 515517, 0.802853 >0.010328/ 

DATA  D1,D2,D3/1. 432788 >0.189269 >0.001308/ 

B=NFl/2. 

A=NF2/2. 

C 

C  CALCULATE  NORMAL  QUANTILE. 

C 

T  =SQRT ( ALOG( 1 . /( ALPHA**2 .  ) ) ) 

NUMER=C0+C1*T+C2*T**2 

DENUM=1 . +D1*T+D2*T**2 . +D3*T**3 . 

XP=T-( NUMER/DENUM ) 

C 

C  CALCULATE  THE  INVERSE  F-DISTRIBUTION . 

C 

LAM=( XP**2-3 . )/6 . 

W2=  (l./(2.*B-l.  ) )-( l./( 2.*A-1 .  )) 

W4=  (l./(2.*A-l.  ) )♦( 1 ./( 2 . *B-1 .  )) 

H=2.*W4**< -1.  ) 

Hl  =  (  LAM-M5./6.  )-( 2./( 3.*H  ) )  ) 

H3=(  XP*(  H+LAM )** . 5  )/H 

H=H3-W2*W1 

FINVER=EXP(2*H) 

WRITEt  3,100  )  NF1,NF2> ALPHA, FINVER 

100  FORMAT (  'O','  NF1= ' ,13 > ' NF2= ' >13 > ' ALPHA= ' , F5 . 4 , * FINVERSE= ' , F6 . 3 ) 

STOP 

END 

The  following  table  shows  some  values  obtained  using  this  approximation, 
together  with  exact  values  of  critical  values  of  the  F  distribution.  It  can  be  seen  in 
Table  8  that  this  approximation  provides  adequate  accuracy  for  the  computation  of  the 
critical  value  of  the  F-tests. 
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TAB  LI:  8 

APPROXIMATION  AND  EXACT  CRITICAL  VALUES  OP  I7- TESTS 


a 

nl 

n2 

approx 

exact 

.05 

4 

4 

6.581 

6.39 

4 

6 

4.601 

4.53 

4 

10 

3.525 

3.48 

10 

20 

2.35 

2.35 

20 

30 

1.932 

1.93 

6 

4 

6.359 

6.16 

10 

6 

4.112 

4.06 

20 

6 

3.944 

3.87 

30 

20 

2.041 

2.04 

30 

30 

1.841 

1.84 

.10 

4 

4 

4.150 

4.11 

4 

6 

3.218 

3.18 

4 

10 

2.648 

2.61 

10 

20 

1.939 

1.94 

20 

30 

1.668 

1.67 

6 

4 

4.046 

4.01 

10 

6 

2.952 

2.94 

20 

6 

2.864 

2.S4 

30 

20 

1.739 

1.74 

30 

30 

1.607 

1.61 
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APPENDIX  B 

THE  INVERSE  CENTRAL  T  APPROXIMATION 


We  can  use  the  asymptotic  expansion  for  the  Inverse  t-CDF,  as 
follows.  Let  Q  be  the  complementary  CDF  of  standard  normal,  and  let  P 
be  the  central  t  CDF  [Ref.  4:  p.  949]. 

If  P(t  |  v  )  =  l-2p  and  Q(x  )  =  p 

then 

tp~  xp  +  (9l  (xp  )/  v)+  (92  (xp  )  /  v2  ) +( 9 3  (xp  )  /  v3  )+(  g4  (xp  )  /  v4  ) 
where  (x)=. 25(x^  +x), 

g2  (x)=( l/96)(5x^  +16x3  +3x), 
g3  (x)=(l/384)(3x7  +19x^  +17x3  -15x),  and 
g4  (x)=( 1/92160) ( 79x9  +776x7  +1482x5  -1920x3  -945x). 

The  program  we  used  in  the  computation  of  the  critical  value  of  the 
t-test  is  shown  below. 


*  * 

*  APPROXIMATION  TO  THE  T-  INVERSE  DISTRIBUTION.  * 

*  CALCULATE  T-INVERSE  GIVEN  DOF  ALPHA.  * 

*  Q(  T-ALPHA  |  DOF  )=  ALPHA.  * 

*  * 


^  *  A  A7\  a  a  AAA  a  a  a  A  A  AAA  a  a  a  a  A  A  A  A  A  A  A  A  A  A  A  A  A  AAA  A  AA  AAAAA  AA  AAaAAA  AAAA  A  7\A  AA  A  A  AA  A  AA  A 

c 

C  set  initial  conditions. 

C 

REAL  ALPHA ,NUMER,DENUM, LAM, XP,G1,G2,G3,X,TINVER 
INTEGER  N 

DATA  N, ALPHA/9,. 025/ 

DATA  CO  ,C1 ,C2/2 .515517,0. 802853,0. 010328/ 

DATA  Dl, 02,03/1.432788, 0.189269, 0.001308/ 

C 

C  CALCULATE  NORMAL  QUANTILE. 

C 

ALPHA= . 05 

T=SQRT ( ALQG( 1 . /<  ALPHA**2 .  ) ) ) 

NUMER=CG+C1*T  +C2*T**2 

DENUM=1 . +D1*T+D2*T**2 . +D3*T**3 . 

XP  =T - ( NUME  R/DE  NUM ) 

C 

C  CALCULATE  THE  INVERSE  T-DISTRIBUTIQN . 

C 

X=XP 

G1=.25*(X**3+X) 

G2  =  ( 1./96.  )*(5.*X**5+16.*X**3+3.*X) 

G3=< 1./384. )*(  3 . *X**7+19 . *X**3+17. *X**3-15 . *X ) 

G4=( 1./92160. )*( 79 . *X**9+776 . *X** 7+1482 . *X**5-1592 . *X**3-945.*X ) 

TP=X+( Gl/N )  +  ( G2/( N*N  ) )  +  ( G3/( N**3 ) )+(G4/(N**4)) 

TINVER=TP 

HRITE( 3,100)  N, ALPHA, TINVER 

100  FORMAT ( '0* , '  D0F=' ,I3,4X, ' ALPHA= ' ,F5.3,6X, 'TINVERSE=' ,F9.5) 

STOP 

END 

The  following  table  shows  some  values  obtained  using  this 
approximation  together  with  exact  values  of  critical  values  of  the  t 
distribution.  It  can  be  seen  in  Table  9  that  this  approximation 
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provides  adequate  accuracy  for  computing  the  critical  regions  of  the 
t-tests. 


TABLE  9 

APPROXIMATION  AND  EXACT  CRITICAL  VALUES  OF  T-TESTS 


DOF 

a  (  Type 

I  error  ) 

.  05 

.  01 

exact 

approx 

exact 

approx 

2 

2.  919 

2.  867 

6.  964 

6.  475 

4 

2.  132 

2.  126 

3.  746 

3.  701 

6 

1.  943 

1.  942 

3.  142 

3.  130 

8 

1.  833 

1.  833 

2.  821 

2.  818 

10 

1.  812 

1.  823 

2.  764 

2.  761 

12 

1.  783 

1.  783 

2.  681 

2.  680 

14 

1.  761 

1.  761 

2.  624 

2.  624 

16 

1.  746 

1.  746 

2.  583 

2.  583 

18 

1.  734 

1.  734 

2.  552 

2.  552 

20 

1.  725 

1.  725 

2.  528 

2.  528 

25 

1.  708 

1.  708 

2.  485 

2.  485 

35 

1.  690 

1.  690 

2.  438 

2.  438 

45 

1.  690 

1.  690 

2.  412 

2.  413 

65 

1.  669 

1.  669 

2.  385 

2.  385 

85 

1.  663 

1.  663 

2.  371 

2.  371 
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APPENDIX  C 

AN  APPROXIMATION  OF  THE  NORMAL  CDF 


Let  Z  be  the  standard  normal  density  function  and  let  P  be  the 
correspondi ng  CDF  [Ref.  4:  p.  939].  Then 

P(X)  =  j  Z(t)dt  ~  l-Z(x)(a2t+a2t^+a3t3)+  £  (x) 

where  t=l/(l+px) 

|  £  ( x) |  <  .000001 

p  =. 33267,  and  a3  =  . 4361836 

a2  =  —  • 1201676 
a3  =  .9372980 
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APPENDIX  D 

THE  COMPUTATION  OF  NONCENTRAL  F  DISTRIBUTION 


***********************  INFORMATION  ********************************** 


*  * 

*  THE  OBJECTIVE  OF  THIS  PROGRAM  IS  TO  CALCULATE  NON-CENTRAL  F  * 

*  DISTRIBUTION  (CDF, F-INVERSE  ).  * 

*  THE  CALCULATION  METHOD  IS  NORMAL  APPROXIMATION.  * 

*  * 

**********************  VARIABLE  DEFINITION  *************************** 
*  * 

*  FVAL  :  X-VALUE  DIMENSION  * 

*  FCDF  :  NON-CENTRAL  F  CDF  DIMENSION.  * 

*  FINVER  :  GIVEN  X  VALUE,  FIND  INVERSE  VALUE.  * 

*  FLAM  :  NON-CENTRAL  PARAMETER  * 

*  FX  :  TO  FIND  INVERSE,  GIVEN  X  VALUE.  * 

*  FSTAR  :  X  VALUE  VARIABLE.  * 

*  NF1  :  DEGREE  OF  FREEDOM  1.  * 

*  NF2  :  DEGREE  OF  FREEDOM  2.  * 

*  CDF  :  NORMAL  DISTRIBUTION  FUNCTION.  * 

*  P  (  Z  <=  X)  * 

*  CONST  :  1  /  SQRT  (  2  *  PI  )  IN  NORMAL  DISTRIBUTION  FUNCTION.  * 

*  B  :  THE  VALUE  OF  'X'  IN  NORMAL  DISTRIBUTION.  * 

*  * 


************************************************************************ 

C 

C  SET  INITIAL  CONDITIONS. 

C 

REAL  FVAL( 100  ) , FCD F ( 100  ) , FINVER , FLAM , FMAX , FX , FINIAL 
DATA  NF1 ,NF2 , FLAM , FMAX/2  >4,75.0,90./ 

DATA  Al, A2, A3, P/. 4361836, -.1201676,. 9372980, .33267/ 

DATA  FSTAR/2S. 128/ 

WRITE( 3,24)  NF1 >NF2 > FLAM 
WRITE (3,21) 

WRITE( 3,22) 

WRITE (3,23) 

CONST  =  1.0  /  SQRT  (2.0  *  3.1415927) 

C 

C  CALCULATE  THE  NORMAL  APPROXIMATION  . 

C 

Cl=( NF1*FSTAR  )/( NF1+FLAM ) 

C2=l-( 2 ./( 9.*NF2 ) ) 

C3=(2.*(NF1+2.*FLAM) )/( 9. *( NF1+FLAM )**2 . ) 

CC1=(C1**( 1./3.  )  )*C2-( 1-C3 ) 

CC2  =  ( ( C3  +  ( 2 . /( 9.*NF2  )  )*( Cl )**( 2 ./3.  )  )**.5) 

FX=CC1/CC2 
CDF  =  0.0 
IF( FX  . EQ . 0  )  THEN 
CDF  =  0.5 

ELSE  IF(FX.GT.O)  THEN 
T=1 ./( 1 . +P*FX ) 

XSQU=( FX*FX  )* . 5 

ZPD  F =CONST*EXP ( -XSQU ) 

CDF=1-ZPDF*( A1*T  +A2*T**2+A3*T**3 ) 

ELSE 

FX=-FX 

T=1 ./( 1 . +P*FX ) 

XSQU=(FX*FX )*.5 
ZPDF=CONST  *EXP ( -XSQU ) 

CDF=ZPDF*( A1*T+A2*T**2+A3*T**3 ) 

END  IF 

FVAL( FSTAR )=FSTAR 
FCDF( FSTAR  )=CDF 

WRITE! 3,250  )  FSTAR ,, FCDF ( FSTAR  ) 

250  FORMAT! '0 ' ,F7.3,8X,F8.5 ) 

21  FORMAT (  'O' ,40(  )) 

22  FORMAT (  ' 0 ' ,2X, 1 "X" -VALUE  1 ,8X, ' F( Z  <=  X)') 

23  FORMAT ( 'O' ,40( '_' )) 

24  FORMAT! ' 0 ' , ' DOFl  = ' ,12 ,3X, ' DOF2= ' ,12 ,3X, 'NON-CENTRAL  PARAMETERS 
&' ,F8.4) 

STOP 

END 
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APPENDIX  E 

THE  COMPUTATION  OF  NONCENTRAL  T  DISTRIBUTION 


INFORMATION  *  *  *  **  **  *  *  *  *  *  *  *  *  ***  **  *  *  **  ******  **  ** 


*.  * 

*  THE  OBJECTIVE  OF  THIS  PROGRAM  IS  TO  CALCULATE  NON-CENTRAL  T  * 

*  DISTRIBUTIONS  (  CDFs )  * 

*  THE  CALCULATION  METHOD  IS  NORMAL  APPROXIMATION.  * 

*  (REF.  SEVERO  &  ZELEN'S  (1960)  )  * 

*  * 

**********************  VARIABLE  DEFINITION  XXXXXXXXXXX XXXXXXXXXXX ^ X  X 
X  X 

*  TVAL  :  X-VALUE  DIMENSION  * 

*  TCDF  :  NON-CENTRAL  T  CDF  DIMENSION.  * 

*  TLAM  :  NON-CENTRALITY  PARAMETER  * 

*  TX  :  TO  FIND  INVERSE,  GIVEN  X  VALUE.  * 

*  TSTAR  :  X  VALUE  VARIABLE.  * 

*  NF  :  DEGREE  OF  FREEDOM  * 

*  CDF  :  NORMAL  DISTRIBUTION  FUNCTION.  * 

*  P  (  Z  <=  X)  * 

*  CONST  :  1  /  SQRT  (  2  *  PI  )  IN  NORMAL  DISTRIBUTION  FUNCTION.  * 

*  B  :  THE  VALUE  OF  'X'  IN  NORMAL  DISTRIBUTION.  * 

*  * 


************************************************************************ 

c 

c  SET  INITIAL  CONDITIONS. 

REAL  TVAL( 1DD),TCDF( 100 ) ,TINVER , TLAM, TMAX,TX,TINIAL 
DATA  NF ,TLAM/24 ,33 . D78/ 

DATA  A1 , A2 ,A3 ,P/ . 4361836 ,- . 1201676 , . 9372980 , . 33267/ 

DATA  TSTAR/24 . 4949/ 

WRITE( 3,24)  NF ,TLAM 
WRITE( 3,21 ) 

WRITE (3,22) 

WRITE (3,23) 

CONST  =  1.0  /  SQRT  (2.0  *  3.1415927) 

CALCULATE  THE  NORMAL  APPROXIMATION  . 

C 1  =  ( 1 . -( 1 . /( 4 . *NF )  )  )*TSTAR-TLAM 
C2=SQRT( 1 . +( TST AR**2/( 2 . *NF  ) )  ) 

TX=C1/C2 

CDF  =  D.O 
IF( TX  . EQ . 0  )  THEN 
CDF  =  0.5 

ELSE  IF(TX.GT.D)  THEN 
T=1 ./( 1 . +P*TX ) 

XSQU= ( TX*TX )* . 5 
ZPDF  =CONST  *EXP ( -XSQU ) 

CDF=1-ZPDF*( A1*T+A2*T**2+A3*T**3  ) 

ELSE 
TX=-TX 

T=1 ./( 1 . +P*TX ) 

XSQU= ( TX*TX )* .  5 
ZPDF=CONST*EXP( -XSQU) 

CDF=ZPDF*( A1*T  +A2*T**2+A3*T**3 ) 

END  IF 

TVALC TSTAR )=TSTAR 
TCDF( TSTAR )=CDF 

WRITE ( 3,250 )  TSTAR , ,TCDF( TSTAR ) 

250  FORMAT ( ' 0 1 ,F7. 3 ,8X,F8.5  ) 

21  FORMAT ( 'O' ,40(  )) 

22  FORMAT ( '0 ' ,2X, ' "X" -VALUE ' ,8X, ' F( Z  <=  X)') 

23  FORMAT (  'O' ,4D(  '_'  ) ) 

24  FORMAT C  'D ■ , 'DOF  = ' ,12, 3X, ’NON-CENTRAL  PARAMETER= 

&' ,F8.4) 

STOP 

END 
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APPENDIX  F 

THE  POWER  OF  ANOVA  PROGRAM  LIST 


* 

*  POWER  ALGORITHM  (  ANALYSIS  OF  VARIANCE  )  * 

*  * 

*  * 

*  DIRECTED  BY  :  PROFESSOR  DONALD . R . BARR  * 

*  * 

*  WRITTEN  BY  :  HUR,  SEONG  PIL  JULY  1986  * 

*  DEPARTMENT  OF  O.R.  NAVAL  POSTGRADUATE  SCHOOL  * 

*  * 

x  * 

*  PROGRAM  USES  THE  POWER  PLOT  TO  TEST  FOR  ANALYSIS  OF  VARIANCE.  * 

*  (A)  ONE-WAY  ANOVA  * 

*  (1)  SAMPLE  SIZE  VS  POWER  * 

*  (2)  NUMBER  OF  TREATMENT  VS  POWER  * 

*  (3)  ALPHA-LEVEL  VS  POWER  * 

*  (4)  NONCENTRALITY  PARAMETER  VS  POWER  * 

*  (B)  MULTI-WAY  ANOVA  * 

*  (1)  NUMBER  OF  REPLICATES  VS  POWER  * 

*  (3)  ALPHA-LEVEL  VS  POWER  * 

*  (4)  NONCENTRALITY  PARAMETER  VS  POWER  * 

*  NOTE  :  subroutine  PLOTT  is  NONIMSL  subroutine  library.  * 

*  * 
*********************************************************************** 

REAL  NN( 200),POWER( 200  ), NON, FINV,NONCEN,FPOW,FX, ALPHA, CDF, CONST 
REAL  MAXDEL ,MINDEL ,INCDEL , STDS, PICON ,MINSD ,MAXSD ,INCSD ,MAXALP 
REAL  MINALP ,INCRAL 

INTEGER  MFI  20 ) ,MAXN, MINN, INCRN ,MAXK , MINK ,INCRK ,JJ >K ,N 
CHARACTER*!  ANS 

C _ 

C 

C  START  ONE-WAY  OR  MULTI-WAY  ANALYSIS  OF  VARIANCE. 

C 

C _ 

500  PRINT*, *DO  YOU  WANT  TO  USE  ONE-WAY  ANOVA  (Y/N)?' 

READ( 5,7)  ANS 
7  FORMAT! Al) 

IF! ( ANS .  EQ . ' N '  ) . OR . ( ANS .  EQ . ' N '  ) )  GO  TO  1000 

PRINT*,' DO  YOU  WANT  TO  PLOT  N(  *  OF  OBSERVATION  PER  TREATMENT) 

&  VS  POWER  !  Y/N)?' 

READ! 5,7)  ANS 

I F( ( ANS . EQ . ' Y ' ) . OR . ( ANS . EQ . 1 Y ' ) )  GO  TO  10 

PRINT*, 'DO  YOU  WANT  TO  PLOT  K  (  NUMBER  OF  TREATMENT  )  VS  POWER 
&  (  Y/N)?' 

READ! 5,7)  ANS 

I F ( ( ANS . EQ . 1 Y '  ) . OR . ( ANS . EQ . 1 Y ' ) )  GO  TO  20 

PRINT*, 'DO  YOU  WANT  TO  PLOT  ALPHA-LEVEL  VS  POWER  (Y/N)?' 

READ! 5,7)  ANS 

IF! ( ANS . EQ . ' Y '  ) . OR . ( ANS . EQ . '  Y '  ) )  GO  TO  30 

PRINT*, 'DO  YOU  WANT  TO  PLOT  NONCENTRALITY  VALUE  VS  POWER  (Y/N)?' 
READ! 5,7)  ANS 

IF! ( ANS.EQ. 'Y ' ) .OR. ( ANS.EQ. 'Y ' ) )  GO  TO  40 
PRINT*, 'DO  YOU  WANT  TO  PLOT  AGAIN  (Y/N)?' 

READ! 5,7 )  ANS 

IFMANS.EQ.  'Y'  ).  OR.  (ANS.EQ.  'Y'  ))  GO  TO  500 
GO  TO  999 

C _ 

C 

C  ANALYZE  NUMBER  OF  SAMPLE  SIZE  PER  TREATMENT  VS  POWER  IN  ONE-WAY 
C  ANALYSIS  OF  VARIANCE. 

C 

C _ : _ 

10  PRINT*, 'NUMBER  OF  K(  #  OF  TREATMENT  )?' 

READ  (5,*)  K 
PRINT*, 'ALPHA-  LEVEL?' 

READ  (5,*)  ALPHA 

PRINT*,'  STANDARDIZED  SQUARED  DISTANCE  VALUE  ?' 

READ  (5,*)  STDS 
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roonooo  oooooo  oooooo  ooooo  noooo 


PRINT*,'  MAXIMUM  N  VALUE  (  THE  MAXIMUM  VALUE  ON  X-AXIS  )?' 

READ  (5,*)  MAXN 

PRINT*,' MINIMUM  N  VALUE  (  THE  MINIMUM  VALUE  ON  X-AXIS  )? 

S  (CONDITION  : N  MUST  BE  MORE  THAN  2  )' 

READ  (5,*)  MINN 

PRINT*, 'INCREMENT  N  VALUE  ?' 

READ  (5,*)  INCRN 
WRITE( 6,95) 

95  FORMAT ( IX , ' DOF1  DOF2  ALPHA-LEVEL  F-INVERSE  NONCENTRAL  POWER') 
NF1=K-1 
JJ  =  1 

DO  110  N=MINN, MAXN, INCRN 
NF2=K*( N-l ) 


COMPUTE  CENTRAL  F-INVERSE  GIVEN  D0F1,D0F2  AND  ALPHA-LEVEL. 


FINV=FINVER  ( NF1 ,NF2, ALPHA ) 


COMPUTE  NONCENTRALITY  PARAMETER  AS  SAMPLE  SIZE  CHANGE. 


NONCEN=N*STDS 


COMPUTE  POWER  USING  NONCENTRAL  F-DISTRIBUTION  (CDF)  GIVEN  DOF1, 
DOF2, CENTRAL  F-INVERSE  AND  NONCENTRALITY  VALUE. 


POW  =  FPO( ALPHA ,NF1 ,NF2 , FINV ,NONCEN  ) 

WRITE (6,15)  NF1,NF2, ALPHA, FINV, NONCEN,POW 
15  FORMAT (  *0' ,I5,3X,I5 ,2X, 2F10 . 5 , IX, F12 . 2 , IX, F10 .5  ) 
NN( JJ  )=N 
POWER( JJ)=POW 
JJ=JJ+1 
110  CONTINUE 
JJ=JJ-1 

PRINT*,' NN  POWER' 

DO  76  N  =  1 , JJ 

WRITE (6,75)  NN( N  )  >POWER( N  ) 

75  FORMAT ('O',  'NN= ' ,F8.3,4X> 'POWER= ' ,F10.5  ) 

76  CONTINUE 


PLOT  NUMBER  OF  SAMPLE  SIZE  PER  TREATMENT  VS  POWER  IN  ONE-WAY 
ANALYSIS  OF  VARIANCE  USING  LIBRARY  SUBROUTINE. 


CALL  PLOTT(NN, POWER, JJ,0) 

PRINT*,' DO  YOU  WANT  TO  PLOT  AGAIN  (Y/N)?' 
READ  (5,7)  ANS 

I F ( ( ANS . EQ . ' Y '  ) . OR . ( ANS . EQ . ' Y '  ) )  GO  TO  500 

GO  T 71  999 


ANALYZE  NUMBER  OF  TREATMENT  VS  POWER  IN  ONE-WAY  ANALYSIS  OF  VARIANCE 


PRINT*, ■#  OF  OBSERVATIONS  PER  TREATMENT  <N=  )  ?' 

READ  (5,*)  N 

PRINT*, 'ALPHA  LEVEL  (REAL  VALUE)  ?' 

READ  (5,*)  ALPHA 

PRINT*, 'WHAT  STANDARDIZED  SQUARED  DISTANCE  VALUE  ?• 

READ  (5,*)  STDS 

PRINT*, 'MAXIMUM  K  (#  OF  TREATMENTS)  RANGE  ?’ 

READ  (5,*)  MAXK 

PRINT*, 'MINIMUM  K  (*  OF  TREATMENTS)  RANGE  (K  MUST  BE  MORE  THAN  2)?' 
READ  (5,*)  MINK 
PRINT*,' INCREMENT  K  VALUE  ?' 

READ  (5,*)  INCRK 
WRITE( 6,195 ) 

195  FORMAT! IX, ' D0F1  D0F2  ALPHA-LEVEL  F-INVERSE  NONCENTRAL  POWER') 
JJ  =  1 

DO  111  K=M INK, MAXK, INCRK 

NF1=K-1 

NF2=K*( N-l ) 

FINV=FINVER  ( NF1 ,NF2 , ALPHA ) 

NONCEN=N*STDS 

POW=FPO( ALPHA, NF1,NF2, FINV, NONCEN) 

WRITE( 6,115)  NF1 ,NF2 , ALPHA , FINV, NONCEN ,POW 
115  FORMAT ( 'O' ,I5,3X,I5,2X,2F10 .5,1X,F12 . 2 , 1X,F10 .5 ) 

NN( JJ  )=K 
POWER( JJ  )=POW 
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JJ=JJ+1 
111  CONTINUE 

PRINT*,' NN  POWER' 

DO  176  N=1,JJ 

WRITE (6,175)  NN( N  ) , POWERC  N ) 

175  FORMAT ('O',  'NN= ' , F8 . 3 ,4X, ' POWER= 1 , F10 . 5 ) 

176  CONTINUE 

CALL  PLOTT (NN, POWER, JJ,0) 

PRINT*, 1  DO  YOU  WANT  TO  PLOT  AGAIN  (Y/N)?' 

READ  (5,7)  ANS 

IF( ( ANS. EQ. '  Y '  )  .OR . ( ANS . EQ. ' Y '  ) )  GO  TO  500 
GO  TO  999 

C _ 

C 

C  ANALYZE  ALPHA-LEVEL  VS  POWER  IN  ONE-WAY  ANALYSIS  OF  VARIANCE. 

C 

C _ 

30  PRINT*, ' *  OF  OBSERVATIONS  PER  TREATMENT  (N=  )  ?' 

READ  (5,*)  N 

PRINT*,'*  OF  TREATMENTS  ( K=  )?' 

READ  (5,*)  K 

PRINT*,'  STANDARDIZED  SQUARED  DISTANCE  VALUE  ?' 

READ  (5,*)  STDS 

PRINT*, 'MAXIMUM  ALPHA  RANGE  ?' 

READ  (5,*)  MAXALP 

PRINT*, 'MINIMUM  ALPHA  RANGE  ?' 

READ  (5,*)  MINALP 

PRINT*, 'INCREMENT  ALPHA  VALUE  ?' 

READ  (5,*)  INCRAL 
WRITE( 6 ,295 ) 

295  FORMAT ( IX , ' D0F1  D0F2  ALPHA-LEVEL  F-INVERSE  NONCENTRAL  POWER') 

JJ  =  1 
NF1=K-1 
NF2=K*( N-l ) 

DO  211  ALPHA=MINALP, MAXALP, INCRAL 
FINV=FINVER  ( NF1 ,NF2 , ALPHA  ) 

NONCEN=N*STDS 

POW  =  FPO( ALPHA ,NF1 ,NF2 ,FINV ,NONCEN  ) 

WRITE( 6,215)  NF1 ,NF2 , ALPHA, FINV ,NONCEN,POW 
215  FORMAT ( 'O' ,15 ,3X,I5 , 2X , 2F10 . 5 , IX, F12 . 2 , IX , F10 . 5 ) 

NN( JJ)=ALPHA 
POWER( JJ )=POW 
JJ=JJ+1 
211  CONTINUE 
JJ=JJ-1 

PRINT*, 'NN  POWER* 

DO  276  I  =1,JJ 

WRITE (6,275)  NN( JJ ) ,POWER( JJ ) 

275  FORMAT ('O',  'NN= ' ,F8.3,4X, ' POWER= ' ,F10.5 ) 

276  CONTINUE 

CALL  PLOTT(NN, POWER, JJ,0) 

PRINT*, 'DO  YOU  WANT  TO  PLOT  AGAIN  (Y/N)?' 

READ  (5,7)  ANS 

IF( ( ANS. EQ. ' Y ' ) .OR. ( ANS. EQ. ' Y ' ) )  GO  TO  500 
GO  TO  999 

C _ 

C 

C  ANALYZE  NONCENTRALITY  PARAMETER  VS  POWER  IN  ONE-WAY  ANALYSIS  OF 

C  VARIANCE 

C 

C _ 

40  PRINT*,'*  OF  OBSERVATIONS  PER  TREATMENT  ( N=  )  ?' 

READ  (5,*)  N 

PRINT*,'*  OF  TREATMENTS(  K=)?' 

READ  (5,*)  K 

PRINT*, 'ALPHA-LEVEL  ?' 

READ  (5,*)  ALPHA 

PRINT*, 'MAXIMUM  STANDARDIZED  SQUARED  DISTANCE  VALUE  RANGE  ?' 

READ  (5,*)  MAXSD 

PRINT*,' MINIMUM  STANDARDIZED  SQUARED  DISTANCE  VALUE  RANGE  ?' 

READ  .(5,*)  MINSD 

PRINT*,' INCREMENT  STANDARDIZED  SQUARED  DISTANCE  VALUE  RANGE  ?' 

READ  (5,*)  INCSD 
WRITE (6,495) 

495  FORMAT ( IX, 'D0F1  D0F2  ALPHA-LEVEL  F-INVERSE  NONCENTRAL  POWER') 

JJ=1 
NF1=K-1 
NF2=K*( N-l ) 

DO  411  STDS=MINSD , MAXSD , INCSD 
FINV=FINVER  ( NF1 ,NF2 , ALPHA  ) 

NONCEN=N*STDS 

POW=FPO( ALPHA, NF1,NF2, FINV ,NONCEN ) 

WRITE( 6,415)  NF1,NF2, ALPHA, FINV, NONCEN,POW 
415  FORMAT( 'O' ,15 ,3X ,15 , 2X, 2F10 . 5 , IX, F12 . 2 , IX, F10 . 5 ) 

NN( JJ )=NONCEN 
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POWER! JJ  )=POW 
JJ=JJ+1 
411  CONTINUE 
JJ=JJ-1 

DO  476  I  =  1,JJ 

WRITE! 6 ,475  )  NN( I ) , POWER! I ) 

475  FORMAT!  'O',  'NN= ’ ,F8 . 3  >4X> ' POWER= ' , F10 . 5 ) 

476  CONTINUE 

CALL  PLOTT(NN, POWER, JJ, 0  ) 

PRINT*,' DO  YOU  WANT  TO  PLOT  AGAIN  (Y/N)?' 

READ  (5,7)  ANS 

I F ( ( ANS .  EQ . ' Y '  ) . OR . ( ANS .  EQ . ' Y '  ) )  GO  TO  500 
GO  TO  999 

C _ 

c 

C  ANALYZE  MULTI-WAY  ANALYSIS  OF  VARIANCE. 

C 

c _ 

1000  PRINT*, ' HOW  MANY  WAYS  ANOVA  CLASSIFICATION?' 

READ! 5,* )  M 
DO  600  1=1, M 

PRINT*, 'FOR', I, '-FACTOR  LEVEL  ,HOW  MANY  LEVELS  ?• 

READ! 5,*  )  MF( I ) 

600  CONTINUE 

C _ 

c 

C  ANALYZE  MODEL  WITH  FACTOR-LEVEL  IN  MULTI-WAY  ANOVA. 

C 

C _ 

PRINT*, 'IS  THERE  ANY  INTERACTION  ?(Y/N)' 

READ! 5,7)  ANS 

IF! (ANS.EQ. 'Y'  )  .OR . ( ANS. EQ. 'Y'  ))  GO  TO  3100 

ISUM=0 

ITOT=l 

DO  800  1=1 ,M 

ISUM=ISUM  +  ( MF ( I  )-l ) 

ITOT=ITOT*MF( I ) 

800  CONTINUE 

IERDF=ITOT-ISUM-l 

C _ 

c 

C  COMPUTE  MINIMUM  NUMBER  OF  REPLICATE  PER  TREATMENT. 

C 

C _ 

MINR=MINIR( ITOT ,ISUM ) 

PRINT*, 'AT  LEAST  MINIMUM  NUMBER  OF  REPLICATE  PER  TREATMENT  =',MINR 
PRINT*, 'WHICH  FACTOR  ANALYSIS  (IF  YOU  WANT  3  RD  FACTOR, PRESS  3)?’ 
READ! 5,*  )  I 
NF1=MF ( I  )-l 

PICON  =  FLOAT( ITOT  )/FLOAT( MF! I )  ) 

GO  TO  2000 

C _ 

C 

C  ANALYZE  MODEL  WITH  FACTOR-LEVEL  AND  2-WAY  INTERACTION  IN 

C  MULTI-WAY  ANOVA. 

C 

C _ 

3100  PRINT*, 'ARE  THERE  ONLY  2-WAY  INTERACTION  AND  FACTOR-LEVEL?! Y/N )  • 
READ! 5,7)  ANS 

IF! (ANS.EQ. 'N*  ). OR. (ANS.EQ. 'N*  ))  GO  TO  2150 

ISUM1=0 

ISUM2=0 

ITOT=l 

DO  850  1=1, M 

ISUM1=ISUM1+! MF! I )-l ) 

ITOT=ITOT*MF ( I ) 

850  CONTINUE 

DO  860  1=1, M-l 
DO  870  J=I+1,M 

ISUM2=ISUM2  +  ( MF ( I )-l )*(  MF! J )-l ) 

870  CONTINUE 
860  CONTINUE 

ISUM=ISUM1+ISUM2 
IERDF=ITOT-ISUM-l 
MINR=MINIR( ITOT ,ISUM ) 

PRINT*, 'AT  LEAST  MINIMUM  NUMBER  OF  REPLICATE  PER  TREATMENT  =',MINR 
PRINT*, 'DO  YOU  WANT  TO  USE  ONLY  FACTOR-LEVEL  (Y/N)?' 

READ! 5,7)  ANS 

IF! (ANS.EQ. 'N'  ). OR. (ANS.EQ. 'N'  ))  GO  TO  880 

PRINT*, 'WHICH  FACTOR  ANALYSIS  (IF  YOU  WANT  3  RD  FACTOR, THEN 
&PRESS  3)?' 

READ! 5 ,*  )  I 
NF1=MF ( I )-l 

PICON=FLOAT( ITOT )/FLOAT(  MF!  I ) ) 

GO  TO  2000 
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880  PRINT*, 'WHICH  INTERACTION  ANALYSIS?(  IF  YOU  WANT  2  AND  3  WAY  INTER 
FACTION, THEN  PRESS  2  3  ).' 

READ( 5 >* )  I,J 

NF1=(  MF( I  )-l )*( MF( J  )-l ) 

PICON=FLOAT( ITOT  )/FLOAT( MF( I  )*MF<  J ) ) 

GO  TO  2000 

C _ 

C 

C  ANALYZE  MODEL  WITH  FACTOR-LEVEL  AND  3-WAY  INTERACTION  IN 

C  MULTI-WAY  ANOVA. 

C 

C _ 

2150  PRINT*,' ARE  THERE  ONLY  3-WAY  INTERACTION  AND  FACTOR-LEVEL  ?(Y/N)’ 
READ( 5,7  )  ANS 

IF<( ANS.EQ.  'N*  ).0R.( ANS.EQ. 'N '  ))  GO  TO  2200 

ISUM1=0 

ISUM3=0 

ITOT =1 

DO  871  1=1, M 

ISUMl=ISUMl+<  MF(  I )-l ) 

ITOT=ITOT*MF(  I ) 

871  CONTINUE 

DO  872  1=1, M-2 

DO  873  J=I+1 ,M-1 
DO  874  K=J+1,M 

ISUM3=ISUM3+ ( MF ( I )-l )*( MF(  J )-l )*(  MF(  K )-l ) 

874  CONTINUE 
873  CONTINUE 

872  CONTINUE 
ISUM=ISUM1+ISUM3 
IERDF=IT0T-ISUM-1 
MINR=MINIR( ITOT , ISUM ) 

PRINT*, 'AT  LEAST  MINIMUM  NUMBER  OF  REPLICATE  PER  TREATMENT  =',MINR 
PRINT*, 'DO  YOU  WANT  TO  ANALYSIS  ONLY  FACTOR-LEVEL  ?(Y/N)' 

READ( 5,7)  ANS 

IF( (ANS.EQ. 'N' ). OR. (ANS.EQ. 'N' ))  GO  TO  875 

PRINT*, 'WHICH  FACTOR  ANALYSIS  (IF  YOU  WANT  3  RD  FACTOR, PRESS  3)?' 
READ( 5 ,*  )  I 
NF1=MF( I )-l 

PICON  =  FLOAT ( ITOT )/FLOAT  (  MF(  I ) ) 

GO  TO  2000 

875  PRINT*, 'DO  YOU  WANT  TO  ANALYSIS  3-WAY  INTERACTION  ?(Y/N)' 

READ( 5,7)  ANS 

IF( (ANS.EQ.  'N*  ). OR. (ANS.EQ.  'N'  ) )  GO  TO  500 

PRINT*, 'WHICH  INTERACTION  ANALYSIS?(  IF  YOU  WANT  2  AND  3  AND  4  INT 
&ERACTION  ,THEN  PRESS  2  3  41.' 

READ( 5 ,* )  I , J ,K 

NF1  =  ( MF ( I  )-l  )*( MF( J )-l )*( MF( K )-l ) 

PICON=FLOAT( ITOT  )/FLOAT( MF( I )*MF(  J )*MF( K ) ) 

GO  TO  2000 

C _ _ 

C 

C  ANALYZE  MODEL  WITH  FACTOR-LEVEL  ,2-WAY  INTERACTION  IN  AND 

C  3-WAY  INTERACTION  IN  MULTI-WAY  ANOVA. 

C 

c _ 

2200  PRINT*, 'ARE  THERE  2-WAY  INTERACTIONS  AND  3-WAY  INTERACTIONS  AND 
SFACTER-LEVEL  (Y/N)?' 

READ( 5,7)  ANS 

IF(  (ANS.EQ. 'N' ). OR.  (ANS.EQ. 'N'  ))  GO  TO  999 

I SUM 1=0 

ISUM2=0 

ISUM3=0 

ITOT  =  1 

DO  950  1=1 >M 

lSUMl=ISUMl+( MF( I  )-l ) 

ITOT=ITOT*MF( I ) 

950  CONTINUE 

DO  960  1=1, M-l 
DO  970  J=I+1,M 

ISUM2=ISUM2  +  ( MF(  I )-l )*( MF( J  )-l ) 

970  CONTINUE 
960  CONTINUE 

DO  965  1=1, M-2 

DO  975  J=I+1 ,M-1 
DO  985  K=J+1,M 

ISUM3=ISUM3  +  (  MF(  I  )-l)*(MF(  J  )-!)*(  MF(K)-l) 

985  CONTINUE 
975  CONTINUE 
965  CONTINUE 

ISUM=ISUM1+ISUM2+ISUM3 
IERDF=IT0T-ISUM-1 
MINR=MINIR( ITOT ,ISUM  ) 

PRINT*, 'AT  LEAST  MINIMUM  NUMBER  OF  REPLICATE  PER  TREATMENT  =’,MINR 
PRINT*,' DO  YOU  WANT  TO  ANALYSIS  ONLY  FACTOR-LEVEL  (Y/N)?' 
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READ (5,7)  ANS 

IF( (ANS.EQ. 'N'  ).0R.( ANS.EQ. 'N'  )  )  GO  TO  990 

PRINT*, 'WHICH  FACTOR  ANALYSIS  (IF  YOU  WANT  3  RD  FACTOR, PRESS  3)?' 
READ(5,*)  I 
NF1=MF( I ) -1 

PICON=FLOAT( ITOT )/FLOAT(  MF(  I ) ) 

GO  TO  2000 

990  PRINT*,' DO  YOU  WANT  TO  ANALYZE  2-WAY  INTERACTION  (Y/N)?' 

READ( 5,7)  ANS 

IF< (ANS.EQ. 'N'  ). OR. (ANS.EQ. 'N' ))  GO  TO  995 

PRINT*, 'WHICH  INTERACTION  ANALYSIS  ?(  IF  YOU  WANT  2  AND  3  WAY  INTE 
FACTION, THEN  PRESS  2  3  ).' 

READ( 5  >* )  I,J 
NF1=(MF(I)-1)*(MF( J)-l) 

PICON=FLOAT ( ITOT )/FLOAT( HF( I )*MF( J )  ) 

GO  TO  2000 

995  PRINT*, ' DO  YOU  WANT  TO  ANALYZE  3-WAY  INTERACTION  (Y/N)?' 

READ( 5,7)  ANS 

IF((ANS.EQ. *N'  ). OR. (ANS.EQ. 'N'  ))  GO  TO  500 

PRINT*, 'WHICH  INTERACTION  ANALYSIS?<  IF  YOU  WANT  2  AND  3  AND  4  INT 
&ERACTION  ,THEN  PRESS  2  3  4  ).' 

READ(5,*)  I,  J,  K 

NF1  =  (MF( I  )-l )*( MF( J )-l )*( MF ( K )-l ) 

PICON=FLOAT< ITOT )/FLOAT( MF( I )*MF(  J )*MF ( K ) ) 

GO  TO  2000 

2000  PRINT*, 'DO  YOU  WANT  TO  PLOT  R(  #  OF  REPLICATE  PER  CELL) 

&VS  POWER  (  Y/N)?' 

RE AD (5,7)  ANS 

IF( (ANS.EQ. 'Y‘  ). OR. (ANS.EQ. 'Y'  ))  GO  TO  4100 

PRINT*, 'DO  YOU  WANT  TO  PLOT  NONCENTRALITY  VS  POWER  (Y/N)?' 

READ( 5,7 )  ANS 

IFUANS.EQ.  'Y'  ).  OR.  (ANS.EQ.  'Y'  ))  GO  TO  4200 

PRINT*, 'DO  YOU  WANT  TO  PLOT  ALPHA-LEVEL  VS  POWER  (Y/N)?' 

READ (5,7)  ANS 

IFUANS.EQ. 'Y'  ). OR. (ANS.EQ. 'Y'  ))  GO  TO  2300 
PRINT*, 'DO  YOU  WANT  TO  PLOT  AGAIN  (Y/N)?' 

READ  (5,7)  ANS 

IF( (ANS.EQ.  'Y'  ). OR. (ANS.EQ.  'Y1  ))  GO  TO  500 
GO  TO  999 

C _ 

C 

C  ANALYZE  NUMBER  OF  SAMPLE  SIZE  PER  TREATMENT  VS  POWER  IN  MULTI-WAY 
C  ANALYSIS  OF  VARIANCE. 

C 

C _ 

4100  PRINT*, 'ALPHA  LEVEL?1 
READ  (5,*)  ALPHA 

PRINT*, 'STANDARDIZED  SQUARED  DISTANCE  VALUE  ?* 

READ  (5,*)  STDS 

PRINT*, 'MAXIMUM  R  VALUE  (  THE  MAXIMUM  VALUE  ON  X-AXIS  )?' 

READ  (5,*)  MAXR 

PRINT*,'  MINIMUM  R  VALUE  (  MUST  BE  MORE  THAN ' ,MINR , '  )?' 

READ  (5,*)  MINR 

PRINT*, 'INCREMENT  R  VALUE  ?' 

READ  (5,*)  INCRR 
WRITE( 6,4095) 

4095  FORMAT ( IX, 'DOF1  DOF2  ALPHA-LEVEL  F-INVERSE  NONCENTRAL  POWER') 
JJ  =  1 

DO  4110  N=MINR, MAXR, INCRR 
NF2=N*ITOT-ISUM-l 
FINV=FINVER  ( NF1 ,NF2 , ALPHA ) 

NONCEN=N*PICON*STDS 

POW=FPO( ALPHA ,NF1 ,NF2 ,FINV,NONCEN ) 

WRITE (6,4015)  NF1 ,NF2 , ALPHA , FINV ,NONCEN,POW 
4015  FORMAT ( '0 ' ,I5,3X,I5,2X,2F10.5,1X,F12.2,1X,F10.5 ) 

NN(  JJ)=N 
POWER( JJ)=POW 
JJ=JJ+1 

4110  CONTINUE 
JJ=JJ-1 

PRINT*, *NN  POWER1 

DO  4076  N=1,JJ 

WRITE ( 6,4075 )  NN( N ) ,POWER( N ) 

4075  FORMAT (  *0 1 ,  'NN= 1 , F8. 3 ,4X , ' POWER= ' , F10 .5 ) 

4076  CONTINUE 

CALL  PLOTT(NN, POWER, JJ,0) 

PRINT* , 1  DO  YOU  WANT  TO  PLOT  AGAIN  (Y/N)?1 
READ  (5,7)  ANS 

I F ( (ANS.EQ. 'Y'  ).OR.( ANS.EQ. ' Y '  ) )  GO  TO  500 
GO  TO  999 

C _ 

C 

C  ANALYZE  NONCENTRALITY  PARAMETER  VS  POWER  IN  MULTI-WAY 
C  ANALYSIS  OF  VARIANCE. 

C 
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c _ 

4200  PRINT* , ' HOW  MANY  NUMNER  OF  REPLICATE  (MUST  BE  MORE  THAN ’ ,MINR, 

&'  )?' 

READ( 5,* )  IRR 
NF2=IRR*ITOT-ISUM-l 
PRINT* , ' ALPHA-LEVE  L  ? ' 

READ  (5,*)  ALPHA 

PRINT*, 'MAXIMUM  STANDARDIZED  SQUARED  DISTANCE  VALUE  RANGE  ?• 

READ  (5,*)  MAXSD 

PRINT*, 'MINIMUM  STANDARDIZED  SQUARED  DISTANCE  VALUE  RANGE  ?' 

READ  (5,*)  MINSD 

PRINT*, 'INCREMENT  STANDARDIZED  SQUARED  DISTANCE  VALUE  RANGE  ?' 

READ  (5,*)  INCSD 
WRITE! 6,4495) 

4495  FORMAT ( IX , ' DOF1  DOF2  ALPHA-LEVEL  F-INVERSE  NONCENTRAL  POWER') 
JJ=1 

DO  4911  STDS=MINSD, MAXSD, INCSD 
FINV=FINVER  ( NF1 ,NF2 , ALPHA  ) 

NONCEN=IRR*PICON*STDS 

POW=FPO< ALPHA, NF1,NF2, FINV, NONCEN) 

WRITE( 6,4915)  NF1 ,NF2 , ALPHA , FINV,NONCEN ,POW 
4915  FORMAT ( ' 0 ' ,15 ,3X ,15 , 2X , 2F10 .5,1X,F12 . 2 ,1X,F10 .5 ) 

NN( JJ  )=NONCEN 
POWER( JJ )=POW 
JJ=JJ+1 
4911  CONTINUE 
JJ=JJ-1 

DO  1476  I  =1,JJ 

WRITE* 6,1475)  NN( I  )  ,POWER( I ) 

1475  FORMAT ('O',  *NN= ' , F8 . 3 ,4X, ' POWER= ' , F10 . 5  ) 

1476  CONTINUE 

CALL  PLOTT(NN, POWER, JJ,0) 

PRINT*,' DO  YOU  WANT  TO  PLOT  AGAIN  (Y/N)?' 

READ  (5,7)  ANS 

IF(  (  ANS.EQ. ' Y '  ) .OR. ( ANS. EQ. *Y '  ) )  GO  TO  500 
GO  TO  999 

C _ 

c 

C  ANALYZE  ALPHA-LEVEL  PARAMETER  VS  POWER  IN  MULTI-WAY 
C  ANALYSIS  OF  VARIANCE. 

C 

C _ 

2300  PRINT*, 1  HOW  MANY  NUMNER  OF  REPLICATE  (MUST  BE  MORE  THAN' ,MINR, 

&’  )?' 

READ( 5,* )  IRR 
NF2=IRR*ITOT-ISUM-l 
READ  (5,*)  K 

PRINT*, 'STANDARDIZED  SQUARED  DISTANCE  VALUE  ?' 

READ  (5,*)  STDS 

PRINT*, 'MAXIMUM  ALPHA  RANGE  ?• 

READ  (5,*)  MAXALP 

PRINT*, 'MINIMUM  ALPHA  RANGE  ?' 

READ  (5,*)  MINALP 

PRINT*, 'INCREMENT  ALPHA  VALUE  ?' 

READ  (5,*)  INCRAL 
WRITE ( 6 ,1295  ) 

1295  FORMAT (IX,’ DOF 1  DOF2  ALPHA-LEVEL  F-INVERSE  NONCENTRAL  POWER') 
JJ=1 
NF1=K-1 
NF2=K*( N-l ) 

DO  1211  ALPHA=MINALP, MAXALP, INCRAL 
FINV=FINVER  ( NF1 ,NF2 , ALPHA  ) 

NONCEN=IRR*PICON*STDS 

POW=FPO( ALPHA ,NF1 ,NF2 , FINV ,NONCEN  ) 

WRITE( 6,1215)  NF1,NF2, ALPHA, FINV, NONCEN,POW 
1215  FORMAT (  ' 0 ’ ,I5,3X ,15 ,2X, 2F10 . 5, IX ,F12 . 2 , IX , F10 .5  ) 

NN( JJ  )=ALPHA 
POWER(JJ)=POW 
JJ=JJ*1 
1211  CONTINUE 
JJ=JJ-1 

PRINT*  >  *NN( JJ ) * ,MINK ,MAXK , INCRK 
DO  1276  I  *1,JJ 

WRITE( 6,1275)  NN( JJ  )  ,POWER( JJ  ) 

1275  FORMAT! ' 0 ' >  'NN= 1 ,F8. 3 ,4X, ' POWER= ' , F10 .5 ) 

1276  CONTINUE 

CALL  PLOTT(NN, POWER, JJ,0) 

PRINT*,' DO  YOU  WANT  TO  PLOT  AGAIN  (Y/N)?1 
READ  (5,7)  ANS 

IF( (ANS.EQ.  'Y'  ). OR. (ANS.EQ.  'Y'  ))  GO  TO  500 
GO  TO  999 


999  STOP 
END 
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ououu  ooooo  *******  ooooo  ooooo 


*  APPROXIMATION  TO  THE  F-  INVERSE  DISTRIBUTION.  * 

*  CALCULATE  F-INVERSE  GIVEN  D0F1 ,D0F2 , ALPHA .  * 

*  Q<  F-ALPHA  | D0F1 , DOF 2  )  =  ALPHA.  * 

*  * 


FUNCTION  FINVER(N1,N2,ALP ) 

REAL  ALP,NUMER,DENUM,LAM,XP,YP,A,B,FIN1 
DATA  CO  ,C1 ,C2/2 . 515517 , 0 . 802853 ,0 .010328/ 
DATA  D1 ,D2,D3/1 .432788,0 . 189269,0 . 001308/ 
B=Nl/2. 

A=N2/2 . 


CALCULATE  NORMAL  QUANTILE. 


T =SQRT ( ALOG( 1 ./(  ALP**2.  ) ) ) 
NUMER=C0+C1*T+C2*T**2 
DENUM=1 . +D1*T+D2*T**2+D3*T**3 
XP=T-( NUMER/DENUM ) 


CALCULATE  THE  INVERSE  F-DISTRIBUTION . 


LAM=( XP**2-3 . >/6 . 

W2=  <1./(2.*B-1.  ))-(l./(2.*A-l. )) 

W4=  ' 1 ./( 2 . *A-1 .  )  )  +  ( 1 . /( 2 . *B-1 .  )) 

H=2.*W4**( -1.  ) 

Wl  =  (  LAM+15./6.  )-(  2./(3.*H)  ) ) 

W3  =  ( XP*( H+LAM )**.5 )/H 
W=(W3-W2*W1 )*2 

IFCW.LT.174)  THEN 
FIN1=EXP(  W) 

ELSE 

FIN1=1000 

ENDIF 

FINVER=FIN1 

RETURN 

END 

* 

GIVEN  DOF1 ,DOF2 ,FINV , NON-CENTRAL  PARAMETER,  * 

THEN  COMPUTE  POWER.  * 

(  I.E.  CALCULATE  NON-CENTRAL  F-DISTRIBUTION  (  C.D.F.)  * 

* 


FUNCTION  FPO( ALP,N1,N2, FIN, NON ) 

REAL  NON 
INTEGER  N1,N2 

DATA  A1 ,A2, A3 ,P/. 4361836 , - .1201676 , . 9372980 , . 33267/ 
CONST  =  1.0  /  SQRT  (2.0  *  3.1415927) 


CALCULATE  THE  NORMAL  APPROXIMATION  . 


Cl=( N1*FIN  )/( Nl+NON  ) 

C2=l-( 2 ./( 9  .*N2  )  ) 

C3=( 2.*(N1+2.*N0N) )/( 9 .*( Nl+NON )**2 ) 
CCl=(Cl**(l./3.  )  )*C2-( 1-C3 ) 

CC2  =  (  (  C3  +  (  2 . /(  9 . *N2  )  )*(  Cl  )**<  2./3.  )  )**.5) 
FX=CC1/CC2 


CALCULATE  THE  NONCENTRAL  F-DISTRIBUTION  USING  NORMAL  CDF. 


CDF  *  0.0 

IF(  FX  . EQ. 0 )  THEN 

CDF  =  0.5 

ELSE  IF( FX.GT.O )  THEN 
T=1 ./( 1 . +  P*FX ) 

XSQU=(FX*FX)*.5 
I F ( XSQU . LT . 1 74 )  THEN 

ZPDF=CONST*EXP ( -XSQU) 

ELSE 

ZPDF=0 

ENDIF 

C  ZPDF=CONST*EXP( -XSQU  ) 

CDF=1-ZPDF*(A1*T+A2*T**2+A3*T**3) 

ELSE 

FX=-FX 

T=1 . /( 1 . +P*FX ) 
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XSQU=( FX*FX)*.5 
I F ( XSQU . LT . 174  )  THEN 

ZPDF=CONST*EXP( -XSQU) 

ELSE 

ZPDF=0 

ENDIF 

CDF=ZPDF*( A1*T+A2*T**2+A3*T**3 ) 

END  IF 
FP0=1-CDF 
RETURN 
END 

**************** XXX **************************************************** 
*  * 

*  GIVEN  ITOT  AND  ISUM  ,  THEN  CALCULATE  MINIMUM  REPLICATE.  * 

*  * 
*********************************************************************** 

FUNCTION  MINIR( ITOT, ISUM) 

INTEGER  IR,  ITOT , ISUM ,ICHECK ,MINIR 
IR=2 

10  ICHECK=I R*I TOT -ISUM -1 

I F I ICHECK . LE . 0  )  THEN 
IR=IR+1 
GO  TO  10 
ELSE 

MINIR=IR 

ENDIF 

RETURN 

END 
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APPENDIX  G 

THE  POWER  OF  T-TEST  PROGRAM  LIST 


*********************************************************************** 


*  * 

*  POWER  ALGORITHM  (  T  TEST  )  * 

*  * 

*  DIRECTED  BY  :  PROFESSOR  DONALD. R.  BARR  * 

*  * 

*  WRITTEN  BY  :  HUR,  SEONG  PIL  JULY  1986  * 

*  DEPARTMENT  OF  O.R.  NAVAL  POSTGRADUATE  SCHOOL  * 

*  * 

*  PROGRAM  IS  USED  THE  PRESENTATION  OF  POWER  ALGORITHM  OF  T-TEST.  * 

*  (1)  ONE-SAMPLE  CASE  TEST.  * 

*  .  ONE-SIDED  TEST  (  HA  :  MU>MU0  )  * 

*  .  ONE-SIDED  TEST  (  HA  :  MU<MU0  )  * 

*  .  TWO-SIDED  TEST  (  HA  :  MU  MUO  )  * 

*  (2)  TWO-SAMPLE  CASE  TEST.  * 

*  .  ONE-SIDED  TEST  (  HA  :  MU>MUO  )  * 

*  .  ONE-SIDED  TEST  (  HA  :  MU<MU0  )  * 

*  .  TWO-SIDED  TEST  (  HA  :  MU  MUO  )  * 

*  NOTE  :  Subroutine  PLOTT  is  NONIMSL  subroutine  library.  * 

*  * 


REAL  NN( 200  )  ,POWER( 200 )  ,DF ,NONCDF ,TCDF1 ,TCDF2 ,TINV1 
REAL  TINV,NONCEN,TX, ALPHA, CDF, CONST  ,MAXMU1 ,MINMU1 ,INCMU1 
REAL  MAXDEL ,MINDEL ,INCDEL ,POW,POWl,POW2 ,MUZERO,STDV ,MU1 
INTEGER  JJ,K,N 
CHARACTER*!  ANS 

C _ 


c 

c 

c 

c 

INPUT  USER  OPTIONS. 

100 

PRINT*,' DO  YOU  WANT  TO  TEST  ONE  SAMPLE 
READ( 5,10 )  ANS 

T-TEST( Y/N )? ' 

10 

FORMAT! Al) 

IF( (ANS.EQ. 'N' ). OR. (ANS.EQ. 'N* ) )  GO  TO 
PRINT*,'  SAMPLE  SIZE  ?■ 

READ( 5 ,* )  N 

PRINT*,'  ALPHA-LEVEL  ?' 

READ( 5,* )  ALPHA 

PRINT*,'  MU-ZERO  VALUE  ?' 

READ( 5 ,*  )  MUZERO 

PRINT*,'  STANDARD  DEVIATION  VALUE  ?' 

160 

c 

READ ( 5 ,* )  STDV 

c 

C  ONE -SAMPLE  ONE-SIDED  TEST,  HA  :  MU  >  MUO. 

C 

C _ 

PRINT*,' DO  YOU  WANT  TO  USE  ONE-SIDED  TEST  (Y/N)?' 

RE AD (5,10)  ANS 

I F ( ( ANS . EQ . ’ N ' ) . OR . ( ANS .  EQ . ' N ' ) )  GO  TO  30 

PRINT* , 1  DO  YOU  WANT  TO  TEST  H0:MU=MU0  VS  HA:MU>MU0  (Y/N)?' 

READ( 5,10)  ANS 

IF( ( ANS.EQ. 'N* ) .OR. ( ANS.EQ. 'N' ) )  GO  TO  40 

PRINT*,'  MAXIMUM  MU1  VALUE  (  THE  MAXIMUM  VALUE  OF  X-AXIS  )?' 
READ  (5,*)  MAXMUl 

PRINT*, '  MINIMUM  MU1  VALUE  (  THE  MINIMUM  VALUE  OF  X-AXIS  )? 

&  (CONDITION  :MU1  MUST  BE  MORE  THAN  MU-ZERO  )' 

READ  (5,*)  MINMU1 

PRINT*,'  LARGE  INCREMENT  N  VALUE  ?' 

READ  (5,*)  INCMU1 
N=N-1 

C_ _ 

C 

C  GENERATE  T-INVERSE  VALUE  (  CRITICAL  REGION  VALUE  ). 

C 

C _ 

TINV=TINVER( N , ALPHA ) 

JJ  =  1 
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DO  50  MU1=MINMU1,MAXMU1,INCMU1 


C _ 

C 

C  COMPUTE  NON-CENTRALITY  VALUE. 

C 

C _ 

DF=N 

NONCEN=( * MU1-MUZERO )*SQRT ( DF+1 )  )/STDV 
C _ 

c 

C  GENERATE  CDF  OF  NONCENTRAL  T  DISTRIBUTION. 

C 

C _ 

TCDF=NONCDF*N>NONCEN,TINV) 

POW=l-TCDF 
NN* JJ )=MU1 
POWER*  JJ)=POW 
JJ=JJ+1 

50  CONTINUE 

JJ=JJ-1 

PRINT*,' NN  POWER' 

DO  60  1=1, JJ 

WRITE*  6 ,70  )  NN( I ), POWER* I  ) 

70  FORMAT*  'O'  ,  'NN=  '  >F8 . 3 ,4X,  1  POWER= '  ,F10 .5  ) 

60  CONTINUE 

CALL  PLOTT*NN, POWER, JJ,0) 

PRINT*,' DO  YOU  WANT  TO  PLOT  AGAIN  * Y/N)?' 

READ  *5,10)  ANS 

IF*  *  ANS. EQ. ' Y '  )  .OR. *  ANS. EQ. ' Y ' )  )  GO  TO  100 
GO  TO  80 

C _ 

c 

C  ONE-SAMPLE  ONE-SIDED  TEST,  HA  :  MU  <  MUO. 

C 

C _ 

40  PRINT*, 'DO  YOU  WANT  TO  TEST  H0:MU=MU0  VS  HA:MU<MU0  *Y/N)?' 
READ*  5,10)  ANS 

IF**  ANS . EQ .  ' N '  ) . OR . *  ANS . EQ .  ' N '  )  )  GO  TO  80 

PRINT*,'  MAXIMUM  MU1  VALUE  *  THE  MAXIMUM  VALUE  OF  X-AXIS  )? ' 
READ  *5,*)  MAXMU1 

PRINT*,'  MINIMUM  MU1  VALUE  (  THE  MINIMUM  VALUE  OF  X-AXIS  )? 

&  *  CONDITION  :MU1  MUST  BE  LESS  THAN  MU-ZERO  )  ' 

READ  *5,*)  MINMU1 

PRINT*,'  INCREMENT  MU1  VALUE  ?' 

READ  *5,*)INCMU1 
N=N-1 

TINV=TINVER*  N, ALPHA ) 

TINV=-TINV 

JJ=1 

DO  90  MU1=MINMU1 ,MAXMU1 , INCMU1 
DF=N 

NONC EN  =  *  *  MU1 -MUZE RO  )*SQRT  *  DF  +  1 ) )/STDV 
TCDF=NONCDF*  N ,NONCEN,TINV ) 

POW=TCDF 
NN*  JJ  )=MU1 
POWER* JJ)=POW 
JJ=JJ+1 
90  CONTINUE 
JJ=JJ-1 

PRINT*, 'NN  POWER* 

DO  110  1  =  1, JJ 

WRITE* 6,120)  NN* I), POWER* I ) 

120  FORMAT*  'O',  ' NN= ' , F8 . 3 ,4X , ' POWER= ' , F10 .5 ) 

110  CONTINUE 

CALL  PLOTTtNN, POWER, JJ,0) 

PRINT*,' DO  YOU  WANT  TO  PLOT  AGAIN  (Y/N)?' 

READ  (5,10)  ANS 

IF* ( ANS.EQ. ' Y '  )  .OR. ( ANS. EQ. ' Y '  )  )  GO  TO  100 
GO  TO  80 

C _ 

c 

C  ONE-SAMPLE  TWO-SIDED  TEST,  HA  :  MU  =/  MUO. 

C 

C _ 

30  PRINT*, 'DO  YOU  WANT  TO  USE  TWO-SIDED  TEST* Y/N)?' 

READ, 5, 10)  ANS 

IF* f ANS.EQ. 'N' ). OR. (ANS.EQ. 'N' )  )  GO  TO  80 

PRINT*, 'MAXIMUM  MU1  VALUE  (  THE  MAXIMUM  VALUE  ON  X-AXIS  )? ' 
READ  *5,*)  MAXMU1 

PRINT*,'  MINIMUM  MU1  VALUE  (  THE  MINIMUM  VALUE  ON  X-AXIS  )? 

&  *  CONDITION  :MU1  MUST  BE  LESS  THAN  MU-ZERO  )' 

READ  *5,*)  MINMU1 

PRINT*,'  INCREMENT  MU1  VALUE  ?' 

READ  *  5 ,*  )INCMU1 
N=N-1 
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ALPHA=ALPHA/2 . 

TINV=TINVER( N , ALPHA ) 

JJ=1 

DO  130  MU1=MINMU1 ,MAXMU1 , INCMU1 
DF=N 

NONCEN=(  (  MU1-MUZERO  )*SQRT  (  DF  +  1 .  )  )/STDV 
TCDFl=NONCDF(N,NONCEN,TINV) 

P0W1=1 . -TCDF1 
TINV1=-TINV 

TCDF2=NONCDF(N,NONCEN,TINVl  ) 

POW2=TCDF2 
POW=POWl+POW2 
NN(  JJ)=MU1 
POWER( JJ ) =POW 
JJ=JJ+1 

130  CONTINUE 
JJ=JJ-1 

PRINT*,' NN  POWER* 

DO  140  1=1, JJ 

WRITE ( 6,150  )  NN( I  )  ,POWER( I  ) 

150  FORMAT ('O',  'NN=' ,F8.3,4X, 'POWER= ' ,F10 .5 ) 

140  CONTINUE 

CALL  PLOTT(NN, POWER, JJ, 0 ) 

PRINT*,* DO  YOU  WANT  TO  PLOT  AGAIN  (Y/N)?' 

READ  (5,10)  ANS 

I F ( ( ANS . EQ . ' Y ' ) . OR . ( ANS .  EQ . 1 Y ' ) )  GO  TO  100 
GO  TO  80 

C_ _ 

C 

C  TWO-SAMPLE  ONE-SIDED  TEST,  HA  :  MUX  >  MUY. 

C 

C _ 

160  PRINT* , ' DO  YOU  WANT  TO  TEST  TWO  SAMPLE  T-TEST( Y/N  )? ' 

RE AD (5,10)  ANS 

I F ( ( ANS . EQ . ' N '  ) .  OR .  ( ANS . EQ . ' N '  ) )  GO  TO  80 
PRINT*,'  SAMPLE  SIZE  FROM  POPULATION  1  ?' 

READ( 5,* )  N 

PRINT*,'  SAMPLE  SIZE  FROM  POPULATION  2  ?' 

READ( 5,*  )  M 

PRINT*,'  ALPHA-LEVEL  ?' 

READ( 5,* )  ALPHA 

PRINT*, 'DO  YOU  WANT  TO  USE  ONE-SIDED  TEST  (Y/N)?' 

READ( 5,10  )  ANS 

I F ( ( ANS . EQ.  1 N '  ).OR. ( ANS.EQ. 'N'  ))  GO  TO  170 

PRINT*, 'DO  YOU  WANT  TO  TEST  HO:  MUX=MUY  VS  HA:  MUX>MUY  (Y/N)?' 
RE AD (5,10)  ANS 

IF( (ANS.EQ. ' N '  ). OR. (ANS.EQ. 'N'  ))  GO  TO  180 
PRINT*,'  MAXIMUM  NONCENTRALITY  VALUE  ?  ' 

READ  (5,*)  MAXDEL 

PRINT*,'  MINIMUM  NONCENTRALITY  VALUE?' 

C  &  (CONDITION  : MU1  MUST  BE  MORE  THAN  MU-ZERO  )? ' 

READ  (5,*)  MINDEL' 

PRINT* , '  INCREMENT  NONCENTRALITY  VALUE?' 

READ  ( 5 ,*  )INCDEL 
Nl=M+N-2 

TINV=TINVER( N1 , ALPHA ) 

JJ  =  1 

DO  190  NONCEN=MINDEL, MAXDEL, INCDEL 
TCDF=NONCDF( N1 ,NONCEN ,TINV  ) 

POW=l-TCDF 
NN( JJ  )=NONCEN 
POWER( JJ )=POW 
JJ=JJ+1 
190  CONTINUE 
JJ=JJ-1 

PRINT*, 'NN  POWER* 

DO  200  1  =  1, JJ 

WRITE( 6 ,210  )  NN( I ) ,POWER( I ) 

210  FORMAT* 'O'  ,  *NN= ' ,F8.3,4X, *POWER= ' ,F10.5) 

200  CONTINUE 

CALL  PLOTT( NN, POWER, JJ,0) 

PRINT*, 'DO  YOU  WANT  TO  PLOT  AGAIN  (Y/N)?' 

READ  (5,10)  ANS 

IFUANS.EQ.  'Y'  ).  OR.  (ANS.EQ.  'Y*  ))  GO  TO  100 
GO  TO  80 

C _ 

C 

C  TWO-SAMPLE  ONE-SIDED  TEST,  HA  :  MUX  <  MUY. 

C 

C _ 

180  PRINT*,' DO  YOU  WANT  TO  TEST  H0:MUX=MUY  VS  HA:MUX<MUY  (Y/N)?' 
READ( 5,10  )  ANS 

IFUANS.EQ. 'N'  ). OR. (ANS.EQ. 'N*  ))  GO  TO  80 
PRINT*,'  MAXIMUM  NONCENTRALITY  VALUE  ?  ' 

READ  (5,*)  MAXDEL 


66 


PRINT*,'  MINIMUM  NONCENTRALITY  VALUE?' 

REAO  (5,*)  MINDEL 

PRINT*,'  INCREMENT  NONCENTRALITY  VALUE?' 

REAO  ( 5,*  JINCOEL 
Nl=M+N-2 

TINV=TINVER( N1 , ALPHA ) 

TINV2=-TINV 
JJ  =  1 

DO  220  NONCEN=MINDEL,MAXDEL ,INCDEL 
TCDF=N0NCDF(N1,N0NCEN,TINV2) 

POW=TCOF 
NN( JJ  )=NONCEN 
POWER( JJ)=POW 
JJ=JJ+1 

220  CONTINUE 
JJ=JJ-1 

PRINT*, 'NN  POWER' 

DO  230  1  =  1, JJ 

WRITE ( 6,240)  NN( I  ),  POWER ( I ) 

240  FORMAT ('O',  *NN= ‘ , F8 . 3 ,4X , ' POWER= ' ,F10 . 5 ) 

230  CONTINUE 

CALL  PLOTT(NN, POWER, JJ, 0 ) 

PRINT*,' DO  YOU  WANT  TO  PLOT  AGAIN  (Y/N)?' 

READ  (5,10)  ANS 

I F  U  ANS . EQ . ' Y ' ) .  OR .  ( ANS . EQ . ' Y ' ) )  GO  TO  100 
GO  TO  80 

C _ 

C 

C  TWO-SAMPLE  TWO-SIOEO  TEST,  HA  :  MUX  =/  MUY . 

C 

C _ 

170  PRINT*,' DO  YOU  WANT  TO  USE  TWO-SIOED  TEST (Y/N)?' 
READ( 5,10)  ANS 

IFUANS.EQ.  'N'  ) .  OR .  (  ANS .  EQ .  'N'  ))  GO  TO  80 
PRINT*,'  MAXIMUM  NONCENTRALITY  VALUE  ?  ' 

READ  (5,*)  MAXDEL 

PRINT*,'  MINIMUM  NONCENTRALITY  VALUE?' 

REAO  (5,*)  MINDEL 

PRINT* , ' INCREMENT  NONCENTRALITY  VALUE? ' 

READ  ( 5,*  JINCOEL 
Nl=M+N-2 
A LPHA1= ALPHA/2. 

TINV=TINVER( N1 ,ALPHA1 ) 

JJ  =  1 

00  250  NONCEN=MINOEL ,MAX0EL ,INC0EL 
TCDF1=N0NCDF( N1 ,NONCEN,TINV ) 

P0W1=1 . -TCDF1 
TINV1=-TINV 

TCDF2=NONCOF ( N1 ,NONCEN ,TINV1 ) 

POW2=TCDF2 
POW=POWl+POW2 
NN( JJ  )=NONCEN 
POWER( JJ )=POW 
JJ=JJ+1 
250  CONTINUE 
JJ=JJ-1 

PRINT*, 'NN  POWER' 

00  260  1  =  1, JJ 

WRITE( 6,270 )  NN( I ) ,P0WER(  I ) 

270  FORMAT ('O',  'NN= ' ,F8.3,4X, 'POWER3' ,F10.5) 

260  CONTINUE 

CALL  PLOTT ( NN , POWER , JJ ,0 ) 

PRINT*, '00  YOU  WANT  TO  PLOT  AGAIN  (Y/N)?' 

REAO  (5,10)  ANS 

IFUANS.EQ. 'Y'  ). OR. (ANS. EQ. 'Y' ))  GO  TO  100 
GO  TO  80 
80  STOP 
ENO 


*  * 

*  APPROXIMATION  TO  THE  T-  INVERSE  DISTRIBUTION.  * 

*  CALCULATE  T-INVERSE  GIVEN  OOF  ALPHA.  * 

*  Q(  T-ALPHA  I  OOF  )=  ALPHA.  * 

*  * 


FUNCTION  TINVER(N, ALPHA) 

REAL  ALPHA, NUMER,OENUM, LAM, XP ,G1 ,G2 ,G3 ,X,TINVER 
INTEGER  N 

DATA  CO, Cl, C2/2. 515517, 0.802853, 0.010328/ 

OATA  01,02,03/1.432788,0.189269,0.001308/ 

C _ 

C 

C  CALCULATE  NORMAL  QUANTILE. 

C 

C _ 
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T=SQRT ( ALOG( 1 ./( ALPHA**2.  )  ) ) 

NUME*=C0+C1*T+C2*T**2 

DENUM=1 . +D1*T+D2*T**2 . +D3*T**3 . 

XP=T-( NUMER/DENUM ) 

C _ 

C 

C  CALCULATE  THE  INVERSE  T-DISTRIBUTION. 

C 

C _ 

X=XP 

G1=.25*(X**3+X) 

G2  =  ( 1./96.  )*<5.*X**5+16.*X**3+3.*X) 

G3  =  ( 1./384.  )*(  3 . *X**7+19 . *X#*3+17 . *X**3-15 . *X  ) 

G4=( 1 . /92160 .  )* ( 79 . *X**9+776 . *X**7+1482 . *X**5-1592 . *X**3-945 . *X ) 
TP=X+(G1/N)+(G2/(N*N)  )+(G3/(N**3  )  )  +  (G4/( N**4 )  ) 

TINVER=TP 

RETURN 

END 

***********************  INFORMATION  ********************************** 


*  * 

*  THE  OBJECTIVE  OF  THIS  PROGRAM  IS  TO  CALCULATE  NON-CENTRAL  T  * 

*  DISTRIBUTION  (CDF)  * 

*  THE  CALCULATION  METHOD  IS  NORMAL  APPROXIMATION.  * 

*  * 


**********************  VARIABLE  DECLARATION  ************************** 

FUNCTION  NONCDF(NF>TLAM,TSTAR) 

INTEGER  NF 

REAL  TSTAR,TLAM,TX,  NONCDF >CDF 

DATA  A1 , A2  >A3 , P/ . 4361836 , - . 1201676 , . 9372980  > . 33267/ 

CONST  =  1.0  /  SQRT  (2.0  *  3.1415927) 

C _ 

C 

C  CALCULATE  THE  NORMAL  APPROXIMATION  . 

C 

C _ 

Cl=( 1 . -( 1 ./( 4. *NF  )  )  )*TSTAR-TLAM 
C2=SQRT ( 1 . +( TSTAR**2/( 2 .*NF  )  ) ) 

TX=C1/C2 

CDF  =  0.0 

IF( TX  . EQ . 0  )  THEN 

CDF  =  0.5 

ELSE  IF( TX.GT.O  )  THEN 
T  =  1 ./( 1 . +  P*TX  ) 

XSQU=(TX*TX)*.5 
I F ( XSQU. LT . 174  )  THEN 

ZPDF=CONST*EXP( -XSQU) 

ELSE 

ZPDF=0 

ENDIF 

CDF=1-ZPDF*( A1*T+A2*T**2+A3*T**3  ) 

ELSE 

TX=-TX 

T=1 ./( 1 . +P*TX ) 

XSQU=(TX*TX)*.5 
IF ( XSQU. LT. 174)  THEN 

ZPDF=CONST*EXP( -XSQU ) 

ELSE 

ZPDF=0 

ENDIF 

CDF=ZPDF*( A1*T+A2*T**2+A3*T**3 ) 

END  IF 

NONCDF=CDF 

RETURN 

END 
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